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•SUMMARY 


A  probability  theory  is  developed  to  describe  the 
initiation  of  detonation. 

A  new  solid  state  reaction  rate  theory  is  presented. 
Some  preliminary  applications  are  made  to  drop-height 
experiments . 

Calculations  are  made  to  determine  the  theoretical 
dielectric  constant  and  loss  tangent  of  a  dielectric 
medium  containing  impurities. 

Conclusions  and  recommendations  are  presented  at 


the  end  of  each  chaptor 
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In  a  previous  report1  we  considered  a  simplified  model 
of  the  Initiation  process,.  For  brevity  we  will  term  this 
the  ’’single-point"  model .  In  this  report  we  continue  the 
development  of  this  model  . 

Section  I  reviews  briefly  the  basic  ideas  of  the  model. 
In  Section  II  the  method  of  moments  Is  applied  to  the  fund¬ 
amental  differential  equation  of  the  model.  Section  III 
shows  lurw  the  saddle  point  method  (method  of  steepest 
descent)  may  be  used  to  solve  the  moment  equations.  Section 
IV  presents  an  Improved  solution  of  the  differential 
equation.  This  Improved  solution  Is  still  based  on  the 
method  of  moments  but  permits  one  to  neglect  higher  moments 
with  much  less  error.  The  method  of  Section  IV  gives  a 
completely  adequate  solution  of  the  differential  equation 
and  thus  constitutes  a  complete  mathematical  treatment  of 
the  model.  Section  V  shows  how  the  model  may  be  made  more 
realistic  while  still  remaining  a  single-point  model.  The 
methods  of  Section  V  have  been  used  in  the  actual  numerical 


1.  First  Quarterly  Technical  Report,  Contract  DA-I^-009- 
ENG-J4.I58,  9  October  1959  to  5  February  I960,  prepared 
by  P.E.C.  Corporation. 
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calculations,.  Section  VI  discusses  the  important  question 
of  the  actual  probability  of  initiation  of  explosion  as 
predicted  by  the  "single-point"  model.  Section  VII 
presents  the  results  of  some  numerical  calculations « 
Section  VIII  presents  some  concluding  remarks  and  makes 
recommendations  for  future  Investigation. 


Section  I  .  The  "Single '-Point"  Model  of  Initiation.  Si,*4.#o*0°$ 


L  kJ  U  U  ■  •  v;  A\/JI 


The  basic  Idea  of  the  "si ogle-point”  ’model  is  to 
assume  that  conditions  are  uniform  over  the  very  small 
region  of  the  explosive  In  which  initiation  occurs <> 

Because  reactions  are  occurring,  the  small  region  is  at 
a  higher  temperature,  T,  than  the  temperature  T  of  the 
surroundings.  However,  as  was  discus sod  in  greater 
detail  in  the  earlier  report,  rather  than  assign  a 
definite  temperature  T  to  tho  small  heated  region ,  vre 
Introduce  a  probability  function  P(T,t)  for  the  region. 
Thus  h(T ,t)  dT  is  tho  probability  that  at  time  t  the  temp- 
eratur©  of  tho  region  lies  between  T  and  T+dT.  Very 
briefly  stated,  the  reason  for  this  approach  is  that  when 
very  small  numbers  of  molecules  are  Involved  wo  cannot 
properly  speak  of  the  number  of  reactions  which  occur  in 
a  time  interval  but  can  only  ntato  the  probability  that  a 
reaction  occurs  in  said  time  interval.  Thus  at  any  time 
we  can  never  really  say  how  many  reactions  have  occurred 
and  therefore  xva  cannot  know'  hew  much  the  temperature  of 
tha  region  has  been  raised.  Rather,  we  can  only  state 
the  probability  of  a  certain  number  of  reactions  and  there¬ 
fore  the  probability  that  the  temperature  has  a  certain 


*0# 


v  alue „ 

The  small  region  (the  "single-point"  since  conditions 
are  assumed t  0  be  uniform  over  the  region)  will  lose  heat 
to  the  surroundings  by  thermal  conduction  at  a  rate  given 
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y* < 

fcv  K(T -V..  ;  .  '•'lit-  cnzwsfctvrfc  K  can  be  related  to  the  tberms.1  , 

.  %^o 

canditc  ti  m:  :.y  ;  k.  m*  the  region  (see  earlier  report)  with  'fa'* 

"*  'c 

the  result  r 


?nM?/  3,4 


% 


with  tj  ~  number  of  molecules  in  the  email  region 
jt  radius  of  the  region 
k  -  thorn  a. L  conductivity  of  the  medium 


N  and  ?my ,  of  course  ,  bo  related  by  introducing  the 
volume  occupied  by  a  single  molecule „  In  Section  V  we 
will  inti- educe  n  different  expression  for  K.  This  will  not, 
however ,  iuvG .1:1  unto  any  of  the  dcivolupment  which  follows 
in  Soot  lone  1 ,  11,  UJ  and  XV  „ 


Reaction.-?  arc  assume d  to  occur  in  the  small  region 
with  n  probability  per  unit  time  given  by 


••=  Me' 


/lf:T 


E  is  the  activation  energy  and  A  is  the  frequency  factor 
of  the  Arrhenius  equation,,  H  is  the  number  of  molecules 
in  the  region,,  k  is  Bolt smarm i s  constant 0  It  will  be 
convenient  to  measure  temperature  in  energy  units,  thus 
in  all  that  follows  T  stands  for  k  times  ^tho  temperature „  If 
wo  chons©  E  itse'-f  as  an  energy  unit  and  1st  U  =  NA  then 
we  have  the  simpler  expression 


If  oach  reaction  releases  an  amount  of  energy  Q. 


we  will  assume  this  to  be  immediately  distributed 
uniformly  over  the  entire  region  of  the  "single-point;".. 
This  is  clearly  a  rather  unrealistic  assumption  but  it 
Is  necessary  In  order  to  preserve  the  "single -point 1 
nature  of  the  model 0  This  energy,  Q,  will  cause  T  to 
increase  by  an  amount  !tq"»  "q"  ia  clearly  related  to 

Q,  N  and  C,  the  specific  heat  of  the  medium,.  The  exact 
expression  Is  given  In  Section  V0 

We  must  now  derive  the  equation  satis if led  by 
P(T,t)0  That  is,  we  must  show  how  to  get  P(T,t+dt) 
if  wo  know  P(T,t)„  The  following  derivation  is  exactly^ 
that  presented  in  the  earlier  report  but  is  developed 
here  in  greater  detail „  What  Is  the  probability  that  at 
time  t+dt  the  temperature  lies  between  T  and  T+AT?  This 
probabili  y  can  cloarly  arise  in  two  ways,, 

1)  The  temperature  at  time  t  may  lie  between  T-q 
and  (T-q)+ATo  If  now  a  reaction  occurs  In  the 
time  dt,  the  temperature  will  be  increased  by 
an  amount  q  and  thus  at  t+dt  the  temperature  will 
lie- between  T  and  T+AT„  The  probability  of  this 
event  Is  clearly  given  by 

t?(T~q,t)ATHp(T-q.)dt] 

The  first  bracket  gives  the  probability  of 
temperature  between  T-q  and  T-q+AT  at  the  time  t„ 

e 


The  b.=acJ<©\;  gives  the  prt-babi lity  tial; 

a  rmo i.i ov,  occurs  in  the  time  interval  dt  If 
the  tc.  .\pt  rat  ire  i;.  T~q*  The  product  of  th.;3-:> 
t’.ro  probabilities  p  ivo a  the  probability  of 
e^e  ;1;  [  ].]  : 

2)  The  temperature  at  time  t,  way  be  slightly 
grc  -it «y  than  T  by  just  enough  so  that  in  the 
time  );itervai  dt  r.fc©  heat  conduction  loss  bring.'' 
the  te.r.n/ r- et'i-riA  dot-r.i  to  T,  Lot  TJ  be  the 
tenroor •vb '.i'o  at  tJr.c  i;v  Then  if  the  temperature 
change..-  ‘‘T  boat  cmdoo'i'i nn  alone,,  it  is  clour 
thn ; 

T  “  TJ-K(T’'-*V0)dt 
,/t-  T  ’  =  i  T-KT()dt )  ( l~Edt )  ~1 
~ ( 'j?*»KT{)d t )  ( 1+Kclt } 

”  T+K(T»T0)dt  (1) 

Sow  rio  it- act  on  is  to  occur  in  the  time  interval 
dt  and  t to  probability  of  this  Is  given  by 
I  p'  T’  )dt  n  Sirica  p(T8)  is  Multiplied  by  dt  and 
is  bHo;  a  first  order  small  quantity,  we  may 
replace  its  argument  T4  by  T„  Thus  the  probability 
of  no  reaction  is  given  by  l-p(T) dt  =  This  must 
bo  nultiplleb  by  the  probability  that  the  temp- 
eratu.ro  is  T  :  at  tbs  time  t0  This  is  given  by 
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P(x  I '  AT  Thus  eJ-Jire  (P)  occurs  with  the 
pi-OL'SbX  j  .1  if 

i'l-'pi  1 )  di:  jP(Tl  j  t }  AT  ‘ 

Y/s  :an  row  a  tat  a  that 

F(T,t-s  dt)  AT  =  t  l-p(T)  dt)P{T<  ,t)  AT  '* 

■f  f  P{  T-q  jt )  AT]  [ p( T-q)  dt] 

tie  now  re* it© 

1  { T ,  'IP  o  t )  P  ( T ,  t )  r  ||  dt 

and 

HT',fc)  *  P(  T-Hi{ T“T0 ) dt »t) 

*  P(T,t)  +  K'(T‘  T  )  dt 

dT 

•□.•so  by  (1:1  l  to t ■  o n  1; tut i  ug  the  relation  (1)  between  T®  and 
T  we  9mo  that 

AT  1  =  ( 1+Kdt) AT 

Igoj  hei-it  c undue vio n  produces  a  compression  of 
ueaporatfurc  infcovffnls  by  the  r actor  ~J— _Q  This 
compression  la  easily  seen  In  the  large  time  limit, 

Vf  the  only  changes  occurring  are  due  to  heat  conduction, 
thou  regardless  of  hoi!?  broad  a  distribution,  P(T),  we 
begin  with,  in  a  -very  large  time  it  must  become  a  very 
narrow  peak  at  the  ambient  temperature  T0o 
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If  we  now  make  the  above  substitutions  and  cancel 
the  common  factor  AT,  we  find 

[P+  ££  at]  *  [l-p(T)dt][P+  K(T-T0)dt][l+Kdt] 

3t  0T 

+P(T-q,t)p(T-q) dt 

=  P+dt [ KP+K( T~T  )i!£  -pP] 

°  3T 

•<*  berms  in  (dt)®  of  higher « 


Thus _  _ _ 

H  ®  KF+K<T-T0)||  ~p(T)P{T)ip(T~q)P(T~q) 


(2) 


Except  where  indicated  iri  the  last  term;  the  arguments 
of  P  are  T  and  t#  This  equation  is  the  basic  differential 
equation  of  the  model..  The  solution  of  this  19  considered 
in  the  next  three  sections » 

The  boundary  condition  on  (2)  is  liirjP(T,t)  iS  0  „ 

1--W* 

The  initial  distribution,  P(T,0) ,  must  be  given  and  will 
usually  be  taken  as  fi(T-T0“»q)  where  6  is  the  Dirac  delta 
function o  This  corresponds  to  a  single  reaction  occurring 
at  t  ~  0  thereby  raising  the  temperature  to  T0+q„  Hie 
temperature  of  the  "point”  can  never  drop  below  the 
temperature  T0  of  the  surroundings,,  This  imposes  the 
condition  F(T,t)  =  0  for  T<£T0» 
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Section  !!<,  The  Method  of  Moments  , 

In  order  to  solve  equation  (2)  we  consider  the 
moments  of  P{T)  defined  by 


5* 


,(t)  »  j  TnP(T,t)dT 


jCffl 


■  0 


The  zero  moment,,  M0(t':  =  ^  P(T,t)ciT,  is  the  probability 
that  the  temperature  of  the  "point"  has  some  value  and, 
therefore,  M0(t)  =  1„  The  first  moment  gives  the  expects- 
tion  {or  mean)  value  of  the  temperature  at  the  point.  We 
will  frequently  denote  it  by  0(t)«  Thus 

Hx(t)  F.  e(t)  --  f  T?  ( T ,  t. )  dT  (3) 

v/m 

■lo 


It  will  be  most  convenient  to  redefine  the  moments 
and  measure  them  relative  to  the  mean  temperature,  0„ 
This  is  a  change  from  the  earlier  report.  Wo  define 


M  ( t )  *  ~9)nP(T,t)dT 


o 


(4) 


It  is  clear  that  M0(t)  =  1  but  that  now  Mx{t)  =  0,  The 
mean  square  deviation,  Ms(t) ,  will  be  commonly  denoted 


by  A(t) . 


The  moment  equations  are  now  obtained  by  multiplying 
(2)  by  (T-6)n  and  integrating  from  Tq  to  «*a .  We  will 
treat  in  turn  the  terms  in  equation  (2)„ 


9 


s  .,n  P,p 

-■*-»  (j  1 

m  Bfc 

T„ 


#  *  5  r  ,' 
!  1 

%r 


m  n  J  | 


■nS_  /  rn  a  \  13 

latv  ' 


«o 


dT 


~  f  (T-e)nPdT  +  n~  r  ( T-6) n  XPdT 
dt  JL  dt  J 

"  0 


Tt~  nMn«-l 


d6 

dt 


K  J  (T-9)nPdT  *  KM. 

Jr£  o 


/j***' 

K  /  (T-0)n(T-T  )^“  dT  *  K  f  CT--0)n(T-8+9<-Tn)fidT 

o  aT  Jr,  0  d-L 


o 

K  f  (T-0)nl  l^dT  +  K(0-T  )  f  iT~e)n£LdT 
J  m  dT  °  Jr.  dT 

o  Ao 

<&,■ 
f  1 


K  J  •3~i'(T  -0)rr,'  ,P]  -  P  ( rrKL )  { T*»  8 ) n  1  dT 

11  ST 


Lo 


+  K<0-To)  j  ||T[  (T. e)nP]  -  nPfT-e)"*1!  dT 

■41 c 

=*  -(n+DlO^  -  nK(0-*To)Mn-1 

-K(T0-8)W1P(T0)  -  K(0-To)  (To~0)nP{To) 


'where  we  have  used  the  condition  that  lim  P(T,t)  =  0„ 

T  -At» 


Thus  we  find 


r 


K  \  (T-9) T~T  jS-dT  =  -(oH)DL  -  nK(0~Tn)M  , 
o  3t  n  o  n~l 


Nothing  can  be  done  with  the  term 


<*** 

/< 

nrt 


,  n 


(T-9')  p(T)P{T)dT 


10 


f  ... 

The  term  j  (T-*9)“piT*»q)F(T*-q)  dT  can  be 
To 

letting  T  =  y+q,  It  becomes 


rewritten  by 


T 

I  (y-8+q)np(y )P(y) dy 

40~q 

If  we  remember  that  P ( y )  =  0  for  y-£T  we  can  change  the 
lower  limit  to  TQ  and  If  we  now  replace  y  by  T  we  get 

(  ( T  -  0+q )  np  ( T )  P(  T )  dT 

Jf 

o 

Pitting  all  these  terms  together  we  get  the  moment 
equations. 


dMn 

dT* 


nM. 


rrl- 


ae  _ 


-  -  nKM. 


'dt 


n 


nK(e-T0)M)>1 


+  j  p(T)P(T)[(TM-q-0)n«(T-9)n]dT 

■^o 

The  first  fav/  equations  ares 


n=0 


dM0 


Remomber  that  M0  -  i  hy 
definition. 


n=l  dMx  de  f 

■ST  +  M4  -  *  K(e-T0)M0  +  fp{T)P(T)(3dT 

Jj 


But  if  we  remember  that  1  and  Mj.  =  0  by 

definition,  this  becomes 


||  88  -  K(G-T0)  +  q  f p(T)P(T)dT 

■A0 
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dA 

dt 


2KA 


*  J  p(T 


)F(T)  [qfc+2q(T.»0)  ]  dT 


n~3 


dM3  de 

■—  +  3KM-  +  3A-- 


dt 


*dt 


3K(6-Tn)A 


+  J  p{T)P(T)  [q3+3q*lT«>0)+3q(T~e)®]dT 


The  higher  order  equations  may  be  simplified  by 
5  the  9  eqt 
find  for  n  2 


using  the  9  equation  to  eliminate  ilx.  If  we  do  this  we 

dt 


dM 

—  +  nKMn 
dt  n 


-r. 


nqM^  ,f  l  pPdT  «*  \  pP[  (T-B+q)  n~(T~9) n]  dT 


We  then  have  the  following  set  of  equations  for  the 
moments : 


M  *  1 
o 

Mj.  «  0 

de 


dt 

dt 

dMa 

dt 

etc , 


+  K9  *  KT 


T°  ,  ,J\ 


pPdT 


2KA  *  pP[q*+2q(T«0)]dT 


J 


+  3KMS  »  pP[q3+3q*(T“0)+3q<t-0)s-3qA]dT 


(6) 


We  must  now  show  how  P(T,t)  can  be  obtained  from 
the  moments.  To  do  this  we  develop  P  in  a  series  of 
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Hermits  functions J 


s 

P(T,t)  =  e”s  S  C_{  fc)H„(s)  (7) 

y ?Tta  n~0 

rn  0 

where  s  = 

Hn  —  the  Hermit©  polynomial  of  order  n„ 
a  is  an  arbitrary  parameter  which  will  be  fixed 
In  a  moment. 

Such  a  development  is  always  possible.  Using  Ms 
development,  we  may  calculate  the  moments  of  P„  This 
is  an  easy  matter  if  use  is  made  of  the  orthogonality 
properties  of  the  Hermite  polynomials.  The  results  ares 


Mq  ~  G0  hence  Cfi  -*  :i 

Mi  “  Cx  hence  0X  ~  0  since  =  0 

Ma  s  A  «  a  +  lj.aOs  thus  if  we  put  a  =  A  then  C H  ~  0 , 

Ma  =  6(2a)3/2Ca  «  6(2A)3/2Ca 

or  C3  =  -——73-  m3 
6{2A)3/2 


In  general 


CU  =  — 


(-1? 


m  N. 


with  N 


n  fiT  m  w  •  (n-2m) ! 

_  nr 


n~2m 

2am 


n 


(2A)n/^ 

We  see  that  a  knowledge  of  the  Mss  determines  the  C*s  and 
thus  P(T,t),  We  have 


P(T,t)  « 


[1  + 


6(2A) 3/2 


MaHa{s)  +  (8) 
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T-0 

where  s  -  a  It  la  clear  that  choosing  a  =  a 
r2A 

gives  the  simplest  possible  expression  for  P  in  terms 
of  the  moments,, 

Section  III.  The  Saddle  Point  Method., 

We  now  show  how  the  moment  equations  (6)  may  be 
solved.  The  difficulty  lies  in  handling  the  terms 
involving  p(T)  on  the  right-hand  side  of  the  equations r 
We  note  from  equation  (8)  that  if  M8  is  small  enough  to 
be  neglected,  then  P  assumes  a  very  simple  Gaussian  form 


P(T,t)  = 


!mA 


Thus  in  the  0  equation  we  need  to 


treat  the  following  integral 


#3 


-M*- 

/1mA 


-1/T  -12^)1 


dT 


In  the  A  equation  we  have  the  integral 


r  -  V _  f-  1/T  - 

J '  TSrf  V  °  24 


t  q®+2q( T-9) ] dT 


Since  we  are  neglecting  M3 ,  we  go  no  further  than  the 
A  equation. 

If  A  is  small,  the  integrands  of  both  I  and  J  are 

very  strongly  peaked  at  T  =  9  because  of  the  factor 
_  (T-9)« 

e  2 A  .  It  would  therefore  be  a  good  approximation 
to  treat  the  rest  of  the  Integrand  as  a  constant  equal 

e, 


to  its  value  at  T 

I 


In  this  approximation  I  becomes 
_  au  -i/e  C  -  (T~e)s 

-;sr  6  L  8  V-  « 


We  commit  negligible  error  by  replacing  the  lower  limit 
T0  by  -  and  the  Integral  i s  then  trivial.  Wo  get 


I  “  qDe  "1/0 

In  the  same  fashion 

.  _  Btt  -1/9 
o  -■  q  ue 

These  results  can  be  improved  upon  by  the  saddle 
point  method o  Use  of  the  saddle  point  method  allows  A 
to  be  larger.  Let  us  apply  the  saddle  point  method  to 


I.  We  write 


where 


I  -jflLr-  I  eJ'/A  dT 

f2rcA  / 


P(T)  =  «  A/T  -  i(T-0)s 


If  we  plot  F(T) t  wo  have  roughly  the  following; 
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We  sec  that  A  "n  puukod  .near  T  =  0  and  f  alia  off  very 
rapidly  to  large  negative  values  on  a! that*  side*  Thus 
the  major  contribution  to  the  integral  comas  from  the 
neighborhood  of  the  peak  and  the  smaller  A  is  the  more 
nearly  is  this  true t  The  peak  of  :?  is  at  T  ~  2  and  Z 
may  be  found  by  putting  f'(T)  =  0.  Hence 

?*(T)  -  A/T®  Cf-0) 


and  Z  is  defined  implicitly  by 


or 


A/Z*  (2  3)  h 
Z  =■-  0  A/Za 


We  see  that  the  smaller  A  i;.;,  the  nearer  is  Z  to  0, 

The  saddle  point  method  nmx  consists  of  expanding  F(T) 
about  T  Z« 

F(T)  *«  F(Z)  +  P'(Z)(T~Z)  :•  iP"(Z)(T-*Z)°  +  .... 
If  we  define  K  -  A/Z- 


Gs  =  1  4  2R 

it  is  readily  seen  that 


F(Z)  -  A/Z(l+iR) 

F '  ('/.)  “  0 
F"v  Z)  «  -  G* 

PW(Z)  »  6R/Z 
F(n)(Z)  -  -  C-l)n  £jj„ 
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We  now  write  jf  an 


i  ■=  -7'4-  I 


r 


1/  A  [  F \  Z )  +i  F "{  Z)  {T-  ?.)*••• » „  ] 
e  aT 


T, 


_  qU  F<Z)/A  (  -  Ga(T-7.)®/2A 

~  77^""  e  /  e  d 

K2nA  -A. 


dT 


where  we  have  neglected  higher  terras  In  the  expansion  of 
F(T)  o  Tf  again  wo  replace  the  lower  limit  by  -*  wo 
have  a  simple  Gauss  Integral  and  we  obtain 


1  qi yjlllllh  ,,  CjU  e-  U+4a)/Z 

1  G  G 


We  note  that  for  very  small  A;  Re.rO ,  G  ("^L  and  Z  • 0 
which  gives  us  I  -  qUc "’^/®  os  before. 

Tho  J  integral  In  treated  in  similar  fashion*  We 

wri  te 

u  C  Fi ?) /A 

t: 

’  T* 


A'nA  J,n 


[qa+2q<T-Z+Z-8HT 


U 

tTSTa 


—  [q^q'Z-P)  J  f  eP(T)/A 


dT 


% 


..  -m~ '  ,r 


dT 

!'  2lt  A  «r/rp__ 

The  first  integral  In  the  one  we  have  already  treated 
and  the  second  Integral  vanishes  in  the  approximations 
we  have  made.  Therefore 


j  rS  $2[q+2(Z-0)  ]eF(Z)//A 

L7 
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It  is  not  necessary  to  neglect  the  higher 
powers  of  (T-Z)  in  the  expansion  of  F(T).,  We  present 
below  the  results  to  the  next  order  of  accuracy,,  The 
details  are  rather  involved  and  we  refer  the  reader  to 
Morse  and  Feshbacb  "Methods  of  Mathematical  Physics” 
for  a  general  treatment  of  the  saddle  point  method,. 

In  summary; 


I  *  g2oF(Z)/A[l-..B1+R£!] 


J  '.i  ^-IeF{Z)//At(q+2RZ)U-Ba+Bs)-2(B11+BiS3)] 


G 

where  Z  ■- 
R 

Gfc  - 
Ba 

Ba 

Bi, 


■*  G  4-  L/7? 

“  A/Z  ' 

-  1+2R 

a  ^(3-3R/2) 
}372 


Pi  7)  a  -  A(l+i'ii)/Z 


,  -~-~l ( 1  r,'<- 67 o  GR+?2 .  3RE  -13'r3/8) 


3 


a1* 

R*ZK 

"cP-T" 


B 


3.  G  ” 


13/P-^T~C2 ’6R+RS) 
G1U 


The  B  terms  are  those  which  arise  from  the  higher 
order  saddle  point  treatment  of  the  Integrals „  These 
terms  are  negligible  if  A  is  sufficiently  small. 

It  should  be  mentioned  that  these  same  techniques 
are  still  applicable  if  the  higher  moments  ( Ma ,  etc.) 


(9) 
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are  not  neglected.  It  is  then  necessary  to  use  the 
saddle  point  method  to  a  higher  order  of  accuracy. 

We  are  thus  faced  with  the  problem  of  solving 
equations  (6)  and  (9).,  With  the  neglect  of  M3  thfjse 
become 

+  KO  *  KT  I 
dt  ° 


dA 

dt 


+  2KA  =  J 


(10 


with  I  and  J  defined  by  (9),  Ths  initial  conditions  are 


0(0)  «  T0  +  q 
t\[0)  «  0 


(11 


The  equations  (9),  (10)  and  (1.1)  form  n  closed 
system  which  can  be  integrated  numerically,.  This  system 
was  programmed  for  the  IBM  6r.O  computer  and  twenty  hours 
of  computation  wore  carried  ouc.  These  computations  soon 
indicated  that  although  A  •  0  for  t  *-  0,  it  rather  quickly 
grows  to  such  a  large  value  as  to  invalidate  the  approxi¬ 
mations  made  earlier  on  the  assumption  that  A  was  small,, 
Further,  as  A  grows  in  size  this  is  an  indication  that 
the  higher  moments,  Ma  etc,.,  are  also  becoming  important. 
We  see,  therefore,  that  the  above  techniques,  which  in 
principle  will  yield  a  solution,  are  in  practice  much  too 
tedious  to  be  useful.  In  the  next  section  wo  sec  how  this 
difficulty  can  be  circumvented..  The  techniques  we  have 
developed  In  Sections  .11  and  III  are  stl.U  applicable  but 
•'  r  <  1  *  gh  t  1  ■/  a  1 1  e re  d  f  o 
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o  ^  -t-  4  TV  rl  In  i:  1  nvi 

l’Ov/  o  -i  w  i  i  j-  »  i  -*  •  ■  f.'  ~i  ■ '  -  '--  *--  w  ^  --  *-  -- 

The  equation  wo  wish  to  solve  is 

-  KP  =  K(T-TD)5j  -p(T)PtT)  +  p(T-q)P(T-q)  1  ?.) 

3^  are 

with  the  initial  conditions  P{T?0)  =  $ (T-T0~a) . 

Now  if  there  were  no  reactions  occurring ,  the  initial  £ 
function  distribution  would  mov«j  to  lower  temperature  with  time 
but  ttfould  remain  a  £  function.  In  fact,  if  we  define 
Yit)  to  be  the  solution  of  the  equations 

mr.  v.-T  1  v(n\  -  T  +  <,  mos 

dtj  AU  ’  J  xo'  ■  ~  r  o  i  >-*-  / 

it  is  clear  that  Y(t)  represents  the  temperature  our  "point" 
would  have  if  it  simply  loses  boat  by  conduction  no 
reactions. 

It  is  novr  evident  that  P(Trt)  =  ^~(T“Y(t))  should 
satisfy  (2)  with  the  reaction  terms  omitted?  namely 

||  -  KP  =  K('ivr0)||  (13) 

To  verify  this  note: 

■If  -  -  §  /«-«,  §  -  S'lT-Y) 

Substituting  In(i3)  we  find 

“  H  -  K/(T~Y)  ••=  K(M0)  £\t-Y) 

2C 


9 


or 


[K(T- 


V 


-f 


Vi  .1  -> 

Tt' 


(T-Y)  +  IC  (T  ’Y)  =  0 


Now  using  (12) 

[K(T-To)  *  K(Y-T0)]  S'  ( T -Y )  +  Kf(T-Y)  =  0 

K(T-Y)  Jf'VihY)  +  Kf(T»Y)  «  0 

TJslng  the  property  of  the  JT  function  that 
xjT(x)  -  -  Jf (x)  vro  have 

-  KJ(T-Y)  +  Kj(T-Y)  *  0  qed  > 

Let  us  now  consider  the  offeci,  of  reactions,,  Tho 
JT  function  will  not  only  sink  to  lower  temperature 
because  of  heat  conduction,  it  will  also  decrease  in 
amplitude  due  to  the  occurrence  of  reactions o  This  we 
may  take  care  of  by  writing 

A(t)  £(T-Y(t))  with  A(0)  =  1„ 

As  time  elapses  the  probability  of  a  reaction  Increases 
and  this  means  the  function  P(T)  will  begin  to  grow  a 
little  bump  in  the  neighborhood  of  the  temperature  f0+2q0 
The  probability  at  T0+2q  will  also  be  undergoing  reactions 
thereby  building  up  a  probability  bump  in  the  neighborhood 
of  T0+3q„  Thus  after  a  short  time  we  would  expect  P(T} 
to  have  the  indicated  character. 
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P(T) 


A 


- 


Jffcfci,  at  Y(t) 

j,./  small  'bump  at  «sTf(  tJ  +  q^To+^q 

*/  ^-gmall  bump  atsrY(t)+2q 


This  picture  sxiggests  trying  a  solution  of  the  forms 

P(T,t)  ~  A(t)  ^(T»Y(t))  +  Z  A  (t)  u^CT.t)  (111) 

m~l 

The  Upj^a  will  be  taken  as  normalized  so  that  J^u^dT  -  1. 
Also  ^  is  localised  In  the  neighborhood  T  ^"Y  +  nq. 

Thus  the  problem  becomes  one  of  determining  the  A's  ns 
functions  of  time  and  the  u^s  as  functions  of  T  end  'Vc 
We  now  substitute  (lit)  into  (2)„  Because  P(T) 
is  different  from  zero  only  in  well  separated  regions, 
we  may  consider  Individually  the  various  regions,.  Consider 
first  the  region  In  the  vicinity  of  T  -  Y(t).  What  terms 
In  Eq<,  (2.)  contribute  In  this  region?  It  la  easy  to  see 
that  we  haves 

U-KA  +  p(T)A]  J~(T-Y)  -  [Y  +  K(T-T0)3  Aj(T-Y) 

and  these  terms  should  equal  zero,, 

Now  using  (12)  we  have 


Remembering  that  x  j  (x)  -  -  #  (x)  we  get 

[A  -  KA  +  p{T)A]  &(  T-Y)  -  KA(T-Y)  /(T-Y)  =  0 

T-Y) 

He  nc  a  [A  -  KA  +  p(T)A  +  KA]  S' (T-Y)  *  0 

or 

I A  *  p(T)  A)  S"  ( T-Y)  *  0 

Because  of  the  function  we  may  also  write  this  as 

(A  +  p(Y)A)  f(T-Y)  0 
For  this  to  vanish  we  find 

A  =  -p(Y)A  !  US) 

This  result  could  clearly  have  been  written  down 
Immediately .  Tfc  simply  sayn  that  the  amplitude  of  the  § 
function at  Y{t)  changes  due  to  the  reactions  which  occur® 
Consider  now  the  terms  which  contribute  In  the 
neighborhood  of  Y(t)  +  q„  Those  toms  ares 

A^Ui  +  Ai  Ui  "  KAiUj,  ~  K(  T“=T_ )  Ai^ysl’  “~p(T)  A^u^ 

8  (16) 

+  p(T-q)  Af(T-q-Y) 

Note  that  the  term  p(T-q)  P(T“q)  in  (2)  furnishes  the 
connection  between  the  two  regions,. 


For  m  .  2  we  have  In  general 


V%  *  ‘'tan,,  -  KA^,  -  K»I0HilS  -  pITUj^ 


How  shall  we  solve  these  equations?  Ea.  Clf?)  is  of  course 
no  trouble,  but  how  do  v;e  treat  'CL 6:  and  i!7)?  The  moment 
technique  is  applicable.. 

Let  ua  define:  8m  '  i  T  «ffldT 

and  M  -  f  T- 6  :  u  dT 
nm  *  m  r a 


Consider  now  E<  ,  ■  Wo  x/i]  1  suppress  the  m  i  subscript 

and  write  fc» ,  u  rather  than  MJ;),  0^  ,  u ■_ Multiply  (16) 
by  (T-B)c  and  integrating  we  meet  the  following  terms 

Ai  (  (T  -0  )nudT  »  AXM 

Ax  J’(T-e)®£ulT  Ax  |»  j\?.-e)rjudT  +  j  (T  >e)n=1m0  wdT 


A,  F  f  nGF  .... 

1  I  11  ?  TpSj,',, 

C;  ‘  '  .  ,‘C 


-  Kax  J  ;t  0)  %dT  ••••  -  kaj;)mti 

KAX  |  (T-  •T0)(T-e)nfflJ‘  dT  ‘i  KAj.  j (T  .G+0-To)  ( T-B) dT 

r  *  ■  *t 

f  .  *  • 

*  KAX  j  (T-e)*+1|^  dT  +  (0-To)'’  J(T-8)*1|$  dT  j 

***  "  *i 

=  j"{T-6)nudT  «  n(e-T o)j  (T-eJ^ndT 

r 


KAi  (n+l)Mn  +  nfS-T^M^ 


2k 


At  \  -G) *  T}lfdT  vl  ieh  cannot  bo  s  iv-' pi  i  lied 

*  J 

+  {  ( T  )fi  p.T-q)  J'T -q-Y)<3T  -•  Aiq+Y-e;73  PvY> 

Putting  these  all  toga  the- r  we  find; 

AiMn  Aj  +  nw.j  +  nWn,xe  -  nk(9-T0)AxMml 


<  18) 


.  A-  ((T--0>n  piT)  u:  T)  dT  +  Ap(Y)  f.  Y*q~0)n 


By  definition 


o 


1  Mn  -  0,  Consider  00)  for  various  n» 


n  -  0  A  a  Ap  { Y )  «■•  Ay  {  p'  T)u{T)  JT- 

* 

n^Jt  Ax0  ~  -  K\e-T(;)Ai  Ap Y)  ‘  Y+'u. -9)  -  Ax  j{7 -0)p<T)u(T)f  dT 


n  '•  ?  Denote  MP  by  A 


e' 


A y  A  ■<-  A,  f  A  *  ,?K  a  f'  !-  Ap  ( Y )  Ytq-0) 

i  V  -  'W 


■«  j(M)  PtT)u(T)dT 


We  may  use  the  n  --  0  squat' on  to  olteinate  A  from  the  n»2 
equation. 


Ax  i  A  •>  ?KA  >~  Ap { Y )  j  Y-i-cp  O)^  -  A  j" 

1.  -  "Y  D  '  ) 


Ax  (  ./(T-O)^  •  A‘f  pi  T)-sa{T)  dT 


If  we  assume  the  higher  momenta  ai^e  negligible,  life  are 
led  to  the  following  system  of  equations*. 


2? 


A  “  ■  p  (  V )  A 


Ap(Y)  -  Aj.  (  p  ( T }  li{  T )  dT 


(19a) 


(19b) 


■  K(e»T0)  -  J<T~0)p(T)u‘fT)dT  +  T£  plYHY+q-0)  '7:9c) 

-  2KA  •  j  f(T-e)2  -  pf  T)u(  T)  dT  p(Y)  ./(Y+q-G)2  -  A 

c. 

(I9d) 


Kobe  the  interpretation  of|j.9b)c  Ax  measures  the 

total  probability  of  the  first  bump  {.since  u  la  normalised  to  unit 

probability).  How  does  A  3  change?  It  gains  from  reactions 

occurring  at  the  ,j,  function  {term  Ap(Y)  ]  and  it  loses  because 

reactions  oeoui'  at  the  first  bump  causing  ’'transitions’’  to 

tho  second  bump  [term  -  Ax  |  p(T)'u(T)  f.Tl  „  u  i3  Df  course 

3 till  an  unknown  rime 1 1, on „  Bot/ever  ,  it  can  be  exp andod 

in  Hermito  functions  just  an  P  was  in  Section  XT, 

..  (T^-0)  ,  . 


n(T,t)  a  -  ff 
tr2-n& 


:  ( t. )  H  (%* 
nv ■  1  nnl  jit 


The  coefficients  are  related  to  the  moments It  is 
easy  to  show  that 


Co  "  i 


Gj.  ~  Mx  -  0 

)  Note:  Gj  and  Ca  vanish 

_  Mg-A  „  ly  because  vo  defined  9  as  the 

Gg  =  -jpp  -  0  since  Ms  =  A  (  mean  of  u  and  A  as  the  moan 

J  square  derivation  of  u„ 

ca  = - ~T7o  Ks 

6(2A)3/2 


CL  =  2 


Ns~2K 


/ti  '  K  K?(6~2 K)f  '  o?K 
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Since  .je  assume  Ma  and  higher  arc  n  egligible  wf  get 


u 


(T»-6)a 

2A 


(21) 


The  u  In  (19)  5s  to  be  Interpreted  as  (21)  and  then 
equations  (19)  form  a  closed  system  for  the  quantities 
A,  Ax,  0 1  A,  all  of  which  are  functions  of  t» 

The  saddle  point  method  of  Section  III  may  be 
applied  without  change  to  the  integrals  of  (19) « 
Omitting  the  B  terms  we  have 


r 


j 


p(T)u(T)dT 


IJ  ~(lHR)/Z 
G  6 


j(T-9)p(T)u(T)dT  ~  (Z-0)~  e"(1+*H>/Z 

.  AU  e-(l+iR)/Z 
Z®G 


|  [  (T*'9)sl"A]p(T)  u(T)dT  =  [  Z-e^-A+^gsr] 

U  e-  (1+*R)/Z 
Cr 


(22) 


It  must  be  remembered  that  ws  have  been  looking  at  the 
first  bump  and  that  many  quantities  in  equations  (19) 
and  (22)  need  the  subscript  ”1”. 

The  hi$aer  bumps,  m?2,  are  treated  in  exactly  the 
same  fashion,,  The  only  difference  is  that  the  J  function 
is  no  longer  present  for  m  ^2  and  therefore  those  terms 
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J  rj  (19b),  (19;),  sv/1.5  ( 19d)  which  involve  Y  are 

replaced  by  Integrals, 

We  give  below  the  complete  set  of  equations  for 
the  $  function  and  the  first  two  bumps „  The  equations 
for  all  higher  bumps  are  obtained  from  those  for  m  =  2 
by  simply  changing  subscripts  * 

m=0  A  -  -p(Y)A 

Ax  «  Ap(Y)  -  Ax  (  p(T)Wj.(T)dT 

Ox  “  -k(9A-T0)  -  /T(T  .e1)p{T)ia1(T)dT 

+  l-pCYMY+q-e*) 

rg B» 

Ai  «  -2kAj.  -  j'  [  ( T-'Oj,  Ip  ( T)u  j,  (T) 


m=l 


dT 


% 


A7r>(Y)[(Y^q.J01)8-A1] 


(23) 


The  integrals  are  given  by  (22)  with  9  -  0! , 

Z  -  Zlt  A  A-i,  ff  -  a1?  B  «  R*.  Z,  G,  and  R 
have  the  same  definitions  as  given  in  Section  III,, 


«p 


A8  “  Ax  Jp(T)ui(T)dT  -  As  \  p(T)uffl(T)dT 


,^Tr 


m=2 


98  -  -k(9a-T0)  -  J„(T-0,)p(T)u8(T)dT 


Ax 


To 


<T+q-0a>P(T>«i(T)dT 


-2k  At 


<&'  •‘•o 

f  [  (T-9iJ)2--Aa3p(T)  Uy  (T)dT 

+  I;  )  [(T+c!“0s)8~^lp(T)u1(T)dT 

“  To 

The  integrals  Involving  -a*  are  given  by  (22)  with 
a»  Z  =  Zs,  A  =  As ,  etc.  The  integrals 
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9-9 


involving  Ux  are  slightly  different. 


( 

I  {T+q-0s)p(T)ui{T)dT 
J'*o 


(Zx+q-9a) 


U_ 


o.(i-Hn1)/z1 


h 


t(T+q-ea)*-&8]p(T)TaxtT)dT 

o 


—  +  <Z1+q-02)®-A2] 
2GX® 


e-(l+iB1)/Zx 


The  Initial  conditions  are 

A  --  1 ,  Am  ~  0  for  m  f;- 1. 

~  Ta  +  (m+l)q 

V  ~  0 

Itt 

This  system  in  amenable  to  solution  on  a  computer. 

It  has,  In  fact,  been  programmed  and  preliminary  calcula¬ 
tions  carried  out.  The  results  indicate  that  the  bumps 
remain  •very  narrow  and  therefore  that  the  approximations 
based  on  small  A  are  very  good. 

The  actual  calculations  are,  however,  performed 
with  the  modifications  of  the  next  section, 

Section  V,  The  Expanding  Single  Point  Model, 

As  indicated  in  the  intro due tion >  it  is  possible 
with  minor  modifications  to  convert  the  "single -point” 
model  into  a  much  more  realistic  model.  It  is  the  purpose 
of  this  section  to  indicate  these  changes.  They  are  two 
in  number;  one  serves  to  remove  the  restriction  to  a  fixed 


number  of  molecules  in  the  heated  region  -  thus  we 
term  this  the  "expanding  single  point”  models  The 
other  change  takes  into  account  the  fact  that  when  a 
molecule  reacts  it  is  then  unavailable  for  further 
reactions „ 

If  an  amount  of  energy  Q  is  released  at  a  point  in 
a  medium  at  t  =  0,  then  the  heat  conduction  equation  pre¬ 
dicts  a  spherical  temperature  distribution  given  by 


T  (y  p  t ) 


CL 


+  !'2 


8^S?2C  1  Dt )  3/2  e 

where  G  is  the  heat  capacity  per  unit  volume* 

D  —  ^  k  is  the  thermal  conductivity » 

T0  is  the  initial  temperature  of  the  medium0 
This  temperature  distribution  is  quite  well  approximated 
by  a  distribution  which  is  constant,  from  r  =  0  to  r  — 
l/kDt  and  drops  to  T0  for  r>/)qDto  In  order  to  conserve 
heat  energy  the  constant  value  of  T  must  be 


(2lt) 


T  =  T  + 
o 


jSu 


80  (IpT  ( Dt }  3/2 


Thus  we  approximate  (Sq)  by 


=  To + 


T£.{ 53t 

T?$  iji)t 


i  25) 
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The  heated  region  is  now  seen  to  be  expanding 

! - 

with  &  radius  proportional  to  \?fc,  The  volume  of  and 
therefore  the  number  of  molecules  in  the  heated  region 
grows  pi’oportlonal  to 

The  quantity  K  introduced  in  Section  I  described 
the  change  in  temperature  of  the  ’’point"  due  to  heat 
conduction  and  was  so  chosen  that 


■f  =  -K<ivr0> 


If  we  differentiate  (2E>)  we  find 


dT  s  X  9,  /-3/2D  \ 

dt  c  ( ) 


1  1 
2  T 


-L.  £ 

32at  C 


1 

(Dt) 3/? 

. 


Thus  we  may  identify  the  K  of  the  preceding  four  sections 

with  ~  A  1  ,  in  order  bo  go  to  the  expanding  point  model, 
2  t 

The  fact  that  K  is  now  a  function  of  time  in  no  way 
Invalidates  our  previous  development* 

If  now  a  second  reaction  occurs  wa  again  have  an 
energy  release  of  amount  Q0  We  must  now  assume  that  this 
energy  is  distributed  uniformly  over  a  region  of  radius 
fkDt  in  order  once  again  to  preserve  the  single  point 
character  of  the  model.  It  Is  easily  seen,  from  energy 
considerations,  that  this  will  produce  a  rise  in 
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temperature  in  the  amount 


_ IQ.  1 

Therefore  'the  q  of  the  preceding  sections  is  now  to 
be  taken  as  a  function  of  t  rather  than  a  constant. 

Again,  all  previous  work  is  valid  with  this  change. 

Now  what  changes  are  necessary  to  take  into  account 
the  loss  of  available  molecules  as  reactions  occur.  We 
©an  introduce  a  characteristic  radius,  r0,  such  that 
^itr0  3  =s  vQ,  the  volume  per  molecule.  This  leads  to  a 
characteristic  time  t0  which  is  equal  to  the  time  it 
takes  the  heat  pulse  (25)  to  spread  out  over  the  region 
occupied  by  a  single  molecule.  It  is  now  clear  that  no 
reactions  can  occur  until  t^t  because  only  then  has  the 
heat  pulse  begun  to  move  into  a  region  where  unused 
molecules  are  available.  We  ceo  then  that  the  number  of 
molecules  available  for  reaction  at  any  time  is  propor~ 
tional  to  (t;  |/2-v,/2> .  W©  should  therefore  write  p(T) 
as  Ue-VT  where  TT  =  V(  tV2 ~t03/2)  with  V  ~  A  is 

V 

the  frequency  factor  of  the  Arrhenius  equation. 

These  remarks  are  valid  if  only  one  reaction  has 
occurred.  If,  however,  two  reactions  have  taken  place, 
we  have  used  up  two  molecules  and  we  should  put 
U  as  V(t^/^-2tQ3/^) .  Hence,  in  the  language  of  Section 
IV,  we  see  that  when  we  calculate  with  p(T)  at  the  6 
function  peak,  W8  shovild  set  U  equal  to  V{ t3/2-t03/2) , 
When  we  calculate  with  p(T)  at  the  first  bump,  we 
must  set  U  equal  to  V(  ■t,3/2~?i:n3/2)  . 
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At  the  second  bumf-  U  —  V(  t  V^~3t  3/2)  f  etc«  In  all 
these  expressions  If  t  is  f;uch  as  to  make  U  negative, 
then  D  is  actually  zero. 

We  see  that  the  higher  bumps  are  only  activated 
after  an  appropriate  interval  of  time.  Thus  the  first 
bump  only  begins  to  appear  when  fp- tQ ,  the  second  bump 
when  t  >  2^/-H0,  etc*  Tills  ha3  the  effect  of  changing 
the  initial  conditions  associated  with  equations  (23) » 

It  is  clear  that  the  calculations  with  (23)  begin  at 
t  »  t0<>  For  t<tQ,  A  ~  1  and  —  0  for  all  rn„  "t0n 
is  a  natural  time  unit  and  in  ths  numerical  calculations 
tD  has  been  set  equal  to  one„  The  initial  conditions 
for  the  first  bump  are : 

Ai  ~  0,  &x  “  0,  -  Y{fc0)tq(t0)  when  t  ~  tQ 

The  initial  conditions  for  the  second  bump  ares 

A3  =  0,  0a  =  OittiJ+qCt.,)  ,  «  -MiiL 

2Gxa(ti) 

when  t  =  tj~  2^/^tn 

Similar  condJ.tions  hold  for  all  higher  bumps  „ 

These  changes  seem  to  complicate  the  equations  a 
great  deal  but,  in  actual  fact,  since  the  computations 
must  be  done  on  a  digital  computer  anyway  the  complications 
are  relatively  minor„ 
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Section  VI.,  Probability  of  Initiation. 

We  c obi©  now  to  the  Important  question  of 
determining  the  probability  that  an  initial  reaction 
will  lead  to  initiation  of  explosion.  To  answer  this 
ve  must  consider  the  behavior  of  the  bumps  with  time, 
in  particular1  the  behavior  of  the  bump  amplitude*  » 

If  v/o  look  at  some  bump ,  say  the  mk-J  one  5  we  see  that 
it  is  gaining  probability  from  the  next  lower  (m-1) 
bump  and  losing  probability  to  fcb©  no xt  higher  (m+l) 
bump.,  At  the  same  time  all  bumps  are  sinking  to  lower 
temper  at  ui'o  because  of  heat  conduction.  Kb  on  a  bump  Is 
first  activated  it  will  gain  more  from  below  than  it 
loses  above «  It  will  thus  increase  In  amplitude  for  a 
time.  Eventually  it  will  begin  to  lose  more  above  than 
it  gains  from  below  and  the  amplitude  will  start  to 
decrease,,  Now  the  m--l  bump  eventually  sinks  to  so  low 
a  temperature  that  it  no  longer  feeds  the  bump.  The 
m^h  bump  continues  to  feed  the  m+l  bump  and  now  we  can 
envision  two  possibilities.  The  first  possibility  is 
that  the  mth  bump  continues  to  transfer  probability  to 
the  m+l  bump  so  rapidly  that  its  amplitude,  A^,  goes  to 
zero  before  it  sinks  so  low  in  temperature  that  It  can  no 
longer  undergo  reactions.  The  second  possibility  Is  just 
the  converse,  namely,  the  bump  does  not  empty  itself 

before  Its  temperature  gets  too  low  to  produce  reactions. 


In  this  case  will  approach  some  constant,  non-zero 
value  * 

In  essence,  the  above  paragraph  simply  states  that 

lim  A-  is  either  zero  or  non-zero*  However,  the  above 

remarks  also  make  clear  the  fact  that  if  lim  A_  -  0 

then  lim  Ag  ®  0  for  all  J!>v\a  Also  if  lim  A^  is  non- 
t-ifc®  t-sw» 

zero  then  11m  A®  ia  non- zero  for  all  ss»  Thus  the 

A«s  split  into  two  groups  defined  by  some  integer 

I,  such  that 

lim  A_(  4  0  m  :•&  m 

fcw* 


lim  A  “  0  m  :*  ra 

...... m 


We  see  that  all  bumps  with  m  m  represent  cases  where 
the  chain  reaction  died  out  and  did  not  lead  to  an 
explosion*  Therefore  the  explosion  probability  is 
given  by 


1  -  lim  (A+Ax+Aa+  ttn%)  (26) 

t^s® 

The  procedure  then  is  to  Integrate  equations  (23) 
out  to  the  last  bump  for  which  A  does  not  go  to  zero* 

This  determines  m  and  then  (26)  gives  the  probability 
of  explosion* 
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Section  VII.  Numerical  Calculations 

In  order  to  proceed  with  numerical  calculations,  it 
is  necessary  to  arrive  at  reasonable  estimates  for  the 
parameters  of  the  theory,,  Since  at  this  stage  no  effort 
is  being  made  to  make  a  direct  comparison  with  experiment, 
we  are  only  concerned  that  the  parameters  have  the  proper 
order  of  magnitude. 

As  pointed  out  previously  the  quantity  T,  referred 
to  as  temperature,  is  actually  being  measured  in  energy 
units.  That  is,  T  actually  signifies  the  temperature 
multiplied  by  Boltzmann^s  constant.  The  natural  unit 
for  energy  is  the  activation  ehergy  E  and  is  of  the  order 
of  a  few  electron  volts.  Since  most  explosives  undergo 
thermal  initiation  in  the  range  from  200°  to  ij.00oC  we 


have  taken  T0,  the  ambient  temperature  to  be  of  the 
order  of  600°K.  In  energy  units  this  is  about  l/£o  ev. 
When  divided  by  an  activation  energy  of  about  2  ev.,  this 
gives  a  T0  of  about  1/tpO .  Thus  we  have  calculated  with 
a  T0  of  the  order  0.02. 


The  quantity  q  Is  given  in  Section  V  as 
q  ss  .3.  . . . &  ,  which  we  write  as  W  . 


32n  0  ( Dt }  3/2 


W  . 
3^ 


Thus 


where 


32k  CD- 


Q  Is  the  energy  release  per  reaction 


C  is  the  specific  heat  per  unit  volume 


where  K  is  the  thermal  conductivity. 
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Specific  heats  are  of  the  order  0.5  cal/cm3°C  and 
thermal  conductivities  range  from  10“%  to  10“^ 
cal  cm  sec°C.  If  we  take  K  about  5  x  10'*^  we  find  for 
D  a  value  of  10  "3  cm^/sec.  How  the  characteristic  time 
t  is  defined  as  the  time  required  for  the  heat  pulse 
to  spread  out  over  the  volume  occupied  by  one  molecule. 
According  to  equation  (25)  this  means  /  liDtQ  equals 
a  molecular  radius.  The  unit  cell  for  the  azles  has 
dimensions  of  the  order  6  A0.  Thus  if  we  put  D  =  10~3 
and  =  3  x  10~®  we  find  tQ  2  x  lO"1^  sec.  Thus 

our  time  unit  will  be  of  the  order  of  10”3-3  sec. 

To  find  a  value  for  W  we  proceed  as  follows.  If 
we  measure  time  in  units  of  t  then  W  is  clearly  the 
temperature  rise  produced  by  the  reaction  after  a  time 
tQ,  I.e.,  after  the  energy  of  reaction  has  spread  out 
over  one  molecule.  Thus  W  ~  . & . . —  .  Now 

•2|S(Dto)V20 

3§S<*.)  is  just  the  volume  occupied  by  a  single  molecule 
and  is  of  the  order  10 ~22  cm^.  Q ,  the  energy  released 
per  reaction,  is  of  the  order  10  ev.  and  the  specific 
heat  per  unit  volume  is  about  1.5  x  10^  ev./cm3°C.  These 
give  for  W  about  5000°C.  Converting  to  electron  volts 
we  get  about  .5  ev.  We  now  must  divide  by  the  activation 
energy  E<*^2  ev.  and  we  finally  arrive  at  W‘^'0.25. 

We  have  written  p(T)  as  V{t^/^-t03/2)e*“^/^. 

If  we  measure  time  In  units  of  t  this  becomes 

o 
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V(t-^2-l)e“k/'I(,  Now  when  p(T)  is  written  in  the  form 
••e/rt 

p  =  NAe  ,  as  is  customary  in  the  chemical  literature, 

A  ranges  from  10^^  to  10 ^  see”3'.  We  must  convert  to 
our  time  unit  which  will  change  A  to  range  from  1  to 
10^.  N  is  the  number  of  molecules  involved,  which  in 
our  case  should  be  taken  as  one.  Thus  we  see  that  V 
should  range  from  1  to  lo\  We  have  calculated  with 
V  =  10s  and  103. 

The  following  graphs  present  the  results  for 
calculations  with  W  =  ,23 ,  V  either  10®  or  10^  and 
T0  either  .019  or  .02.  These  calculations  woro  carried 
out  to  the  third  bump  with  the  hope  that  the  value  m 
defined  in  Section  VI  would  be  reached.  This  did  not 
prove  to  be  the  case  as  the  results  show.  The  calculations 
did  not  extend  beyond  the  third  bump  because  the  program, 
as  written,  used  the  full  capacity  of  the  IBM  650  computer. 
However,  the  results  show  that  the  bumps  remain  much 
narrower  than  expected.  That  is,  the  quantity  Aj^, 
which  measures  the  width  of  the  rc^th  bump,  remains  very 
small  throughout  the  calculation.  This  means  that  to 
very  good  accuracy  we  may  ignore  the  width  of  the 
bumps  thereby  greatly  simplifying  the  calculations.  We 
have  almost  completed  a  new  program  which  can  compute 
an  unlimited  number  of  bumps. 

Although  the  calculations  presented  do  not  reach 
the  value  m,  they  do  otherwise  exhibit  exactly  the 
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3.0  "  '  ~h  i/o  u mi+5  o-V  "to  4-To 


the 


■behavior  predicted  in  Section  VI.  The  bump  amplitudes, 

V  start  at  zero,  increase  to  a  maximum  and  then 
decrease  finally  to  a  constant  value. 

Section  VIII.  Conclusions  and  Recommendations. 

The  calculations  prssented  are  extremely  encouraging. 
We  feel  that  the  nearly  completed  program,  designed  to 
calculate  an  arbitrary  number  of  bumps,  will  allow  us  to 
reach  the  m  bumps.  This  will  permit  computation  of 
initiation  probabilities.  We  strongly  urge  the  continua¬ 
tion  of  the  computational  program. 

We  have  Just  recently  made  some  progress  in  an  attempt 
to  obtain  a  completely  analytical  solution  to  the  ’’single- 
point”  model.  Such  a  solution  would  be  an  extremely 
valuable  addition  to  the  numerical  results  as  it  would 
allow  the  general  nature  of  the  "single-point"  model  to 
be  more  adequately  explored.  We  feel  this  to  be  a  very 
promising  area  of  investigation. 

Perhaps  a  few  remarks  concerning  the  significance 
of  the  model  are  in  order.  The  model  is  admittedly 
much  simplified.  We  feel,  however,  that  it  is  a  signifi¬ 
cant  exploration  of  the  probabilistic  aspect  of  the 
initiation  process.  After  the  changes  of  Section  V,  we 
feel  that  the  only  seriously  unrealistic  aspect  of  the 
model  which  remains  is  the  matter  of  assuming  that  the 
energy  Q  of  a  reaction  is  spread  uniformly  over  the 


entire  heated  region.  It  is  clear  that  a  more  realistic 
treatment  of  this  point  would  make  the  probabilistic 
approach  even  more  necessary. 

It  is  our  feeling  that  a  completely  adequate 
treatment  of  the  initiation  process  might  contain  the 
treatment  of  this  report  as  a  bridge  between  a  more 
detailed  treatment  of  the  earliest  stages  of  initiation 
and  the  latter  stages  when  the  heated  region  has  grown 
to  a  macroscopic  hot  spot  for  which  a  probability  approach 
is  unnecessary.  We  have  begun  to  look  at  possible 
approaches  to  the  earliest  stages  of  the  initiation 
process  and  have  some  indication  of  progress  in  this 
direction. 


A  MOISC'JLAR  THEORY  OP 
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Section  I.  Probability  of  Reaction  as  a  Function  of 
Interraolecular  Distance* 

The  potential  energy  may  be  plotted  as  a  function  of 
Interraoleetilar  distance  as  shown  schematically  in  Figure 
1.  In  the  diagram 
Is  the  equilibrium  distance 
and  the  curve  I  gives  the 
potential  energy  as  a 
function  of  r„  If  the  temp¬ 
erature  is  not  high, 
molecules  In  a  solid  will 
vibrate  about  an  equilibrium 
position  in  an  approximately 
simple  harmonic  way  as  will 
be  discussed  in  the  next 
section.  If  the  amplitude  of  vibration  becomes  large 
enough,  the  substance  can  go  Into  new  equilibrium  config- 
uration  with  potential  energy  represented  by  curve  II. 

If  the  amplitude  becomes  larger  there  will  be  a  probability 
of  the  reaction 


N,  +  N, 


3Ne  +  2e 


(1) 


where  electrons  on  the  right-hand  side  are  taken  up  by 
the  metal  Ions . 


The  situation  can  be  illustrated  schematically  by 
means  of  a  series  of  equipotential  energy  contours  as 


shown  in  Figure  2.  In  this  diagram  the  potential  energy  is 
expressed  as  a  function  of  the  -  N8“  separation  and 
also  as  an  average  of  the  N3  -  Na  distances.  The  solid 
equipotential  curves  show  potential  energy  contours  when 
the  2N3  configuration  is  maintained „  The  dotted  equipo¬ 
tential  curves  schematically  represent  the  potential  energy 
contours  whsn  the  nitrogens  become  predominantly  groups, 
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Schematic  of  Equipotential  Contours 
Flguro  2 

Since  the  configurations  have  different  electronic 
structure,  the  two  potential  surfaces  can  co-exist  at 


N 


i,  J 
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each  point  of  this  diagram,,  Tbs  minimum  of  solid-lined 
potential  surface ,  of  course,  corresponds  to  the  equilibrium 
distance  of  Na  ions,  while  the  dotted-line  potential  surface 
does  not  have  a  minimum  except  at  infinity  since  N8  molecules 
are  known  to  have  no  bound  state  (a  molecule).  The  cross- 
section  of  this  diagram  along  the  broken  straight  lines  is 
our  Figure  1,  where  the  curve  II  is  for  the  3^s  configuration 

Reaction  (1)  can  occur  if  one  brings  two  N3  ions 
along  the  curve  I  and  crosses  the  point  where  a».(;r0~r) 

Is  aQ,  This  aQ  ia  called  the  critical  amplitude.  It  is 
possible,  however s  that  even  if  the  Ks"  ions  are  put  together 
with  distance  smaller  than  this  critical  amplitude,  they  come 
back  to  the  2NS  configuration  following' the  curve  I,  They 
have  to  jump  into  the  curve  II  in  order  that  reaction  (l) 
take  place. 

The  probability  that  molecules  jump  from  the  curve  I 
to  the  curve  II  can  be  assumed  to  be  proportional  to  the 
time  that  the  2Na  system  spends  In  the  region  over  the 
critical  amplitude  a0.  The  probability  can  actually  be 
a  rather  complicated  function  of  the  Na""  distance.  The 
above  assumption  corresponds  to  the  simplest  case  that  the 
probability  Is  constant  only  if  the  distance  Is  smaller 
than  this  critical  amplitude. 

The  tunnel  effect  can  produce  the  reaction  even  if 
the  amplitude  of  the  oscillation  Is  smaller  than  the 
critical  value,  but  the  probability  of  a  reaction  through 
such  a  mechanism  must  be  negligible. 


Suppose  aQ  Is  large  so  that  classical  mechanics 
gives  a  good  approximation „  If  the  vibration  is  simple 
harmonic,  the  displacement  x  is  given  by 


x  ~  x  sin  at 
o 


Practiori  of  a  period  spent  in  time  dt  is  dt/T.  Thus  ths 
probability  that  the  particle  be  found  in  range 
x^jd'dx  is 


Pdi 


dt  „ 

T  ' 


dx 


(2) 


Tea  ./x®“X*  211  J x*-x“ 

*  o  S'  o 

Let  ^  s'  x )  be  the  reaction  rate  when  the  particle 
is  at  x.  If  a  reaction  occurred  only  for  x  ap  then 
(x)6(x-a0)  dx  -  S’(a)  »  for  example. 

If  there  is  a  certain  probability  Pdx  of  particle 
being  found  between  x,  and  x*dx  then  Pdx  gives  the 
number  of  reactions  per  unit  time  while  the  particle  is 
in  the  range  x^x+dx»  Thus  the  effective  reaction  rate 
for  the  simple  harmonic  oscillator  is 

P  f'*- 


Ki 


<f(x)Pdx  *  i 
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(3) 


Suppose  the  reaction  rate  is  such  that 


5(x)  ~  <£ : 


if  x  >  a 


if  x  <  a 


Ip 6 


*  0 


(4) 


where  i,  is  a  constant,  me  effective  reaction  rate  K  is, 
from  equation  {3) 


Case  1 „  Model  for  reaction  rate  ${x\ 

Figure  3 

f  is  just  the  fraction  of  n  capable  of  reacting. 


if  xo  <  a0 
if  xo  >  a0 


(6) 
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Case  1,  K3.(n)  as  a  function  of  amplitude  x0» 


Figure  i|. 

F0  In  soma  cases  one  may  expect  the  reaction  rata  to 
be  a  delta  function  as  shown  in  Figure  b » 
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Case  -h  j?(x)  a  delta  function  at  x  “  a0„ 


That  i s 


Figure  5 
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a  3  s  introduced  so  that  may  be  of  the  same  dimension 
o 

as  f'.  In  this  case  we  denote  the  effective  reaction  rate 

"  a. 


hri 


constant  by  Ks. 


;>ne  ooi  aina  T‘..’o  < 


equation  (’3)  again 


!‘/8a06U~a0idt  «  f 


xc  f* *"** 


Integrating  thin  ,  one  has  for  the  effective  reaction 
rate  Ka 


Kj,  ": 

0 

If  x  <  a 

O  0 

Ke  « 

x 

r*. 

§  %  _  £> 

if  x0>a0 

V  -o'ao 

qua  1.3  tatlve 

plot  of  K2 

&m  a  function 

( 8) 


xn  Is  dhown  in  Figure  6f. 


Case  Kj}  as  a  function  of  the  amplitude  x  „ 

Figure  6 

3»  T>ere  can  be  a  case  such  that  the  reaction  probability 
is  a  linear  function  of  tin  velocity  and  occurs  when 
•x.  r;  «,»  so  teat  it,  ir  given  by  the  operator 


'Six)  ‘rsT5(x' 


aQ)v 


,v? 


(9) 


,e '  e  '?  is  ;1  e  period  of  oscillation  and  v  the  velocity 
nA  x„ 


'lie  rate  constant  is  then 


K3  *  ^Tp'°P6U~aoWdx  «  f,,T  f  o)  /^? 

!  I  i/t  _£  .  .,,2 


dx 


•x. 


■Xr 


If  /  {  X^,“XS 


r  r;i  03 
*  3  1  ~» 

•tr 


>  VV 


7'h  (,•  s.. 


E;>  «  0 


i.t‘  x  <  a 
o  o 


K»  "  2$3  it'  K  >  & 


do) 


!'■'.>  this  case  effective  reaction  rate  is  simply  a 
constant  whi-i  the  amplitude  exceeds  a  „  Plot  of  K3  as 
a  function  of  x0,  then  is  siirri.ly  0.3  shown  in  Figure  7. 


Case  v.-  K4i  as  a  function  of.  the  amplitude  x0„ 

Figure  7 

Spotton  II »  Sinst&in  Approximation* 

In  order  to  obtain  the  overall  reaction  rate,  that 
is.,  the  total  number  of  reactions  per  second  per  unit 
mass  of  a  sol  .d  which  is  in  a  thermal  equilibrium,  it  is 


noK  neceaarry  to  calculatt  distribution  of  amplitudes 
x0  at  a  given  temperature., 

In  the  Einstein  approximation  we  regard  the  crystal 
as  consisting  of  If  Independent  oscillators,  all  of  which 
are  vibrat:  r-g  with  the  common  circular  frequency  w0  If 
we  assume  that  the  oscillators  are  harmonic  the  energy 
Is 

S  *•  .W*x£  (11) 

which  is,  in  quantum  mechanics,  nllw  where  n  is  the 
quant  um  number  and  may  take  on  positive  Integer  values <, 


'Inns 


x. 


i  £iM. 

“mo> 


(12) 


The  probability  of  finding  an  oscillator  in  the 
quantum  rjfcf. tj  n  is  given  by 

-r-fMo/kT 
■i  „„  Imo/tcT 


(13) 


Eir«t.-  In  wt  3  able  to  account  for  the  essential  behavior 
of  S)>:'  boat  capacity  of  cryata  1.  by  this  model*  By  means 
of  oh  5  3  mod  a  L  one  car.  now  calculate  an  average  reaction 


"ie<  nu  n  )or  of  reactions  which  occur  In  unit  time 
err  ci."  is  -hen  given  by 


1 


o-*nnm/kT 
-- —  — — — dn 

I...,, 


(lit) 
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where  K\  n)  i*  the  ei so.t5.vo  reaction  rate  calculated  in 
the  preview;  section  expressed  as  a  function  of  the 
quantum  nurb nr  n„  By  using  three  expressions  for  K 
obtained  before  we  can  perform  the  above  calculation 
as  follows : 


Casa 


where 


C*  Klln;e-"fi“/kT 


f 

•'0 


•1 1 


1 1- 


“e/Tj 


-sin”1  1° 

f  i  *o 
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raw  k 
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Cons  idol'  the  portion  of  the  integral  in  sin'"-'  —  . 

xo 


The  expehen  w.al  is  a  maximum  a ;  n  --  nD  and  decreases  very 
rapidly  as  o  increases..  Therefore  an  accurate  expansion 
for  sin~  (n0/n)^  i :■  raquirsci  only  for  n^nQ{  s i .n“-<-n  /n  ^  n/2)  „ 
Now  sin  '.t./2-j)  •='  !.'-yE/,?  so  u/2-y  -  8in“'M  l~y*/2) 

Let  z  =  '-y®/ 2-  y  /ifl-z  ' ,  since  a  ^  1,  let  2  *?  1+z 


Then  n/2- ^  sin'"-1- a  ' 

So  that  ain“:L(nD/n)t  15  n/2-ln-n0)/n  «  n/2“(n-n0)/n0 
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So  Integra)  “cr  K3.  is 
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Letting  ft  •-  integral  in  3  become* 
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where  nQ  -  — ,  b  -  ^ 

—bn 

Consider  the  integral,  I  =|  . — dn»  Let  V  «  n-n, 
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Case  3 

Por  Kg  the  integral  is  simply 
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with  the  same  expression  for  T0  as  before. 

In. order  to  compare  the  theoretical  models  with 
experimental  results,  it  Is  necessary  to  know  how  K  varies 

sit 


with  T  experimentally,,  The  principal  variation  of  K  with 
temperature,  of  course,  will  occur  via  the  Arrhenius 
factor  because  of  the  large  values  obtained  for  EQ  (of  the 
order  of  20-ij.O  kilicalories  per  mole) .  The  "activation 
temperature",  that  is,  the  quantity  E0/R  will  then  vary 
between  10,000°  and  20,000°  Kelvin,  At  the  time  of  this 
writing,  we  have  not  had  an  opportunity  to  examine  all  avail¬ 
able  experimental  data  on  the  reaction  rate  constants  as  a 
function  of  temperature  in  order  to  determine  whether  the 
smaller  effect  associated  with  the  f(T)  in  the  reaction 

rate  constant  1?  »  Af (T)e~®o'^T  can  bo  separated  from  the 

-•E  /RT 

effect  of  the  Arrhenius  term  e  '0/  '  (when  R  rather  than 

K  appears  in  the  exponential,  it  is  understood  that  Eq  Is 
then  expressed  in  calories  per  mole).  Careful  study  of 
experimental  data  should  be  carried  out.  A  brief  discussion 
relating  to  some  experimental  results  for  NH^NO^  will  be 
presented  in  a  later  section. 

There  are  experimental  situations  other  than  a  direct 
determination  of  K,  however,  which  may  enable  one  to  say 
something  about  the  applicability  of  f(T)  to  the  explanation 
of  experimental  results.  Reaction  rate  constant  appears 
in  the  general  conservation  of  heat  equation 

C  =/tV*T  +  Qn& 

3t 

where  C  is  the  specific  heat,  /^is  the  density,  /?  is  the 
thermal  conductivity,  Q  the  heat  of  reaction  per  mole, 
n  the  number  of  moles  per  unit  volume,  and  K  the 
calculated  reaction  rate  constant.  To  reduce  the 


complexity  of  this  equation  and  also  to  define  a  simple 
experimental  arrangement,  one  may  Investigate  solutions 
for  either  an  adiabatic  arrangement  or  for  a  steady  state 
condition. 

The  stationary  state  Is  defined  by 

a  ^eT  =  oriK 

Explosion  takes  place  when  this  equation  is  not  satisfied, 
that  is,  when  a  stationary  state  cannot  exist,  Frank - 
KamenetskI  [Acta  Phys.  Chem.  URSS  10,  365  (1939)1 
considered  this  type  of  problem  with 

K  -  A  exp( -E/kT) 

where  A  is  a  temperature  independent  factor.  We  shall 
Investigate  the  effect  of  the  temperature  dependence  of 
A  as  given  by  the  three  models, 

A  very  simple  experimental  arrangement,  which  in 
principal  should  be  subject  to  exact  mathematical  descrip¬ 
tion,  is  one  in  which  a  gample  of  the  explosive  is  placed 
In  a  temperature  controlled  bath  or  environment.  Both  the 
explosive  and  the  bath  start  out  at  the  same  temper ature , 
Because  of  internal  reactions,  however,  the  explosive  will 
generate  heat  Internally  and  its  temperature  will  rise. 

One  may  then  adjust  the  temperature  of  the  surrounding  bath 
so  that  its  temperature  rises  simultaneously.  In  this  way 
no  heat  Is  lost  from  the  explosive.  That  is,  the  term 
is  equal  to  0  and  the  reaction  takes  place  under 
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adiabatic  conditions  The  heat  flow  equation  to  be 
solved  then  is  simply 

cf  !£  -  Qrif 

In  this  case  the  equation  can  be  integrated  and  the 
temperature  of  the  explosive  determined  as  a  function  of 
time  up  to  detonation.  The  f(T)  term  in  1?  will  yield 
somewhat  different  results  depending  on  the  form  of  f(T)0 
One  may  then  attempt  to  compare  these  theoretical  results 
with  experimental  curves  to  determine  the  importance  of 
f(T)  o 

At  the  time  of  this  writing  we  have  not  been  able 
to  find  experimental  data  for  this  type  of  adiabatic 
experiment.  However,  some  theoretical  curves  have  been 
calculated  and  the  results  will  be  presented  in  the 
next  section. 


sdy  "tote  and  '  AcH  p.bati 
fixper  irnerrcs 


*\'30  H,1  r  ■ '.  f,I>  ..-I..  j>lo«  l'.  a  ? . 

Stoar'Y  ?ha'e  Theory 

In  the  Frank -Kamenatski  the ory  explosion  takes 
place  for  those  values  of  r  and  T  for  which  the  following 
equation  is  satisfied; 

<?  =  S»£  g(T)  e'S/kT 

^  A 

where  E  =  a£M»ri/2 

4 

^  is  a  constant  which  depends  on  the  geometrical  shape  of 
the  explosive,  r  is  the  linear  dimension  of  it,  N  is  the 
number  of  reacting  molecules  per  unit  volume,  Q  is  the 
heat  of  reaction  and  $  is  the  thermal  conductivity.  In 
the  Fr ank-Kame ne t ski  case: 


In  our  Case  is 


g(T)  =  M-_ 
kTB 


_  AT0 
-  -hpr 


g(T) 


In  our  Case  2s 


g(T) 


_  Tp 

(I~e-8/T)T3/2 


In  our  Case  3: 

gCT)  = 


!2p„>*V0 

Tt0(l-e"0/,T)T 


The  difference  among  these  four  cases  appears 
when  the  explosive  material  becomes  very  small.  We 
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notice  that  in  the  Frank -Katnenetski  case  there  exists  an 


absolute  critical  radius  such  that  below  that  radius 
explosion  can  never  take  place  even  at  high  temperature*, 
The  critical  radius  is  obtained  by  solving  for  rmj _  in 
the  Frank -Kane net ski  relationship 


v( 


or 

T* 


e  o 


T0/T 


where 


'o  k 

\jhich  can  be  written  as 
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T„  -  ?,  n  *  SSfl® 
fi  k 


„8  _  i.Ts  Jo/* 

V  —  "Z~  Q 


for  the  critical  radius, 


Thus  for  Frank-Kamenetskl  case,  the  absolute  critical 
radius  Is 


rc  = 


/2F 

V  TvKtav 


T*  V  QilAk 

where  e  is  the  exponential,  and  the  corresponding 
temperature  is 


Tc  «  T0/2 

As  can  be  seen  from  the  general  expression  for  re , 
an  appreciable  change  of  the  critical  radius  r  appears 
when  T  becomes  comparable  to  T0*  By  looking  at  a  curve 
we  can  say  that  such  change  of  r  with  T  occurs  at  about 


T  “  T0/2.5.  If  we  take  this  criterion  we  see  from  the 
formula  for  rE  that  the  critical  radius  at  such  tempera¬ 
ture  is  about  1.1  r  ,  which  means  if  the  change  of 
explosion  temperature  is  observable,  the  absolute  critical 
temperature  is  also  observable. 

In  our  cases  one  can  obtain  the  absolute  critical 
radius  in  the  same  way,  and  one  may  verify  that  they 
occur  at 

Tr  *■  2(TQ+9)  for  Case  i 

Tr  +0)  for  Case  2 

'•  .3  o 

T  «  T  +6  for  Case  3 

co 

Since  the  T-  dependence  of  the  critical  radius  is 

determined  primarily  by  the  exponential  function  , 

the  abovo  criterion  T  ”  T0/2„5  in  order  to  have  an 

appreciable  change  of  r  with  temperature  may  be  applied 

In  a 13  of  the  above  cases  also.  Thus  the  radius  at  which 

the  explosion  temperature  changes  are  approximately 

1,9  ru  for  Case  1 

1  „1  r  for  Case  2 

c 

1  -Ip  ?  for  Case  3 

c 

Experimental  results-*-  show  a  change  of  explosion 
temperature  with  the  dimension  of  explosive.  However,  there 
Is  no  definite  result  for  the  absolute  critical  radius. 

Thus  our  Case  1  is  favorable  over  other  cases  including  the 
original  Frank-Kamenetski  case. 

1.  Pi  P.  -Bowden  and  A.  C.  McLaren,  cited  in  p,  30  of 
Bowden  and  Yoffe«s  "Past  Reactions  in  Solids ",19580 
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The  Adiabatic  Heat  Balance  Equation  is 


where  0  ~  heat  capacity/mole 

Q  =  heat  of  decomposition/mole 
For  the  three  functions , % ( x) s  considered  in  Section  I, 


K  will  be  ox  the  form 

K  =  Af(T)e"T°/T 


where 


f  (T) 


Tn 


and  A  Involves  the  constants  which  vary  with  the  particular 
model.  For  the  present  we  shall  regard  A  a.s  determined 
from  experimental  data  for  a  particular  substance. 

That  Is,  no  attempt  will  be  made  to  evaluate  A  from  more 
fundamental  properties  of  the  solid. 

Tile  hoat  balance  equation  becomes 


'I 


dT  = 

Since  f{T)  i3  a  slowly  varying  function  compared 
T  /*J! 

with  e  0/  ,  we  shall,  at  present,  carry  out  the  integration 
in  T  by  treating  f(T)  as  constant.  Thus  one  obtains 


t-tx) 


T|  bToAx  _  Te  T0/T 
f(Tx)T0  8  '  r(T) To  6 


Thus,  for  each  of  the  three  cases  one  obtains  the 


results : 
Case  1 

t-t2.  = 


QA,t0 


6To//T;i  -  T^{1  -e“e/T)eTo//T 

A  _  1 


2(  5tno)  '*^3/2 


Case  2 


QAgT  f 


As 


..  /«n0* 


Case  3 


t-t-i  = 


Aa 

Git 

-T  /T 

When  the  Arrhenius  equation,  “  Ae  0/  Is  used 
for  K,  the  result  is 


t-t  =  ^ 


T§eTo/Ti  „  T2eTo/^T 
When  the  Eyring  form,  «  At,Te""T°/'r  is  used. 


one 


obtains 


t~tx  “ 


a5qt 


Ti8T»/-  -IaV* 
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The  following  calculation  is  an  exploratory  one 
to  determine  the  nature  of  the  temperature  versus  time 
curve  for  the  adiabatic  decomposition  of  a  typical 
explosive,.  In  the  calculation  it  is  assumed  that  an 

p 

initial  reaction  rate  constant  is  known  .  It  will  be 
assumed  that  the  reaction  rate  constant  is  a  function 


only  of  the  temperature  and  not  of  the  time.  One  may 

arbitrarily  take  the  Arrhenius  value  of  A  for  EDNA  which 

2  11  1 

according  to  Cook  is  10  .  We  will  assume  that  this 

value  of  A  is  determined  from  the  experimental  isothermal 


reaction  rate  constant  at  a  temperature  of  l\27°K*  The- 
constant  in  our  easel,  which  we  will  call  a”  is  then 


given  by 


A"  « 


where  T  Is  the  temperature 


at  which  A  is  determined. 


in  this  case  lj.27°K«  One  may  then  assume  an  adiabatic 
reaction  starting  from  some  arbitrary  initial  temperature, 
sQy  373°K,  and  carry  the  calculation  to  it33°K.  Results 
of  such  a  calculation  are  shown  in  Figure  8  for  two 
arbitrarily  selected  values  of  the  Einstein  temperature 
9  =*  200 °K,  0  «  800°K. 

One  may  note  that  there  is  a  20-30  percent  difference 


Jn  the  predicted  time  from  reaction  initiation  to 

( 


2  ■>  "Isothermal  Decomposition  of  Explosives"  by  M.  A. 

Cook  and  M.  Taylor  Abegg,  Industrial  and  Engineering 
Chemistry,  ij.8,  p.  1090,  June  n 
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detonation  using  the  Arrhenius  reaction  rate  constant  and 
using  ouf  Case  1  reaction  rate  constant «  The  adoption 
of  a  reaction  rate  constant  determined  from  an  isothermal 
experiment  to  describe  the  adiabatic  reaction  may  be 
open  to  question  but  the  curves  do  give  one  a  qualitative 
indication  as  to  how  the  temperature  time  predictions  will 
behave  for  the  various  forms  of  the  reaction  rate  constant. 

Using  the  experimental  data  for  log  K(T)  obtained 
by  Cook  (see  Ref,  2)  for  ^ ,  an  attempt  was  made  to 

determine  whether  our  Case  1  would  fit  the  K(T)  curve  more 
accurately  than  the  simple  Arrhenius  reaction  rate  constant. 
It  was  found  that  no  measureable  improvement  was  obtained 
by  introducing  the  factor  t3/^(  l-©**®/^) „  Since  the  data 
used  in  the  calculation  were  taken  directly  from  the 
graphs  of  Ref.  It  is  possible  this  calculation  may 
be  improved  by  referring  to  the  original  data. 

Conclusions  and  Recommendations 

The  reaction  rate  theory  developed  in  this  chapter 
predicts  a  higher  function  of  the  temperature  in  addition 
to  the  usual  Arrhenius  factor  in  describing  solid  state 
reaction  rate  constants.  Predictions  of  the  theory 
therefore  differ  from  those  made  on  the  basis  of  the 
Arrhenius  or  EyTing  theory. 

It  is  recommended  that  further  work  be  carried  out 
in  order  to  develop  improved  models  for  the  function 
defined  in  the  text.  Further  applications  also  should  be 
made  to  reaction  dependent  phenomena. 


Chapter  III 


PRELIMINARY  CALCULATIONS  RELATING  TO 
EXPLOSION  INDUCED  BY  A  PALLING  BODY 

M.  Mizushlma 

A  shock  wave  can  excite  the  lattice  vibrations. 
If  the  maximum  amplitude  is  xQ  the  energy  of  an 
oscillator  per  atom  is  imsMxQ.  Thus  in  a  simplified 
model  which  assumes  that  n  oscillators  are  excited  by 
the  same  amount  gives  the  total  energy  as 


¥  =  §osMx|? 
d  o 


The  number  n  must  be  proportional  to  the  duration  time 
of  the  excitation  divided  by  the  velocity  of  the  sound  wave 
in  the  explosive.  Suppose  such  shock  wave  is  excited  by 
a  body  which  changes  its  momentum  from  p  to  zero  during 
time  At  we  have 


if 


=  P  =  inkx„  =  InM»8x 
**  M  2  0  2  0 


coAt  <<  1 


If  At  is  much  larger  than  l/co,  we  have 

;  ft  =  wit  Mw8xo  if  1 

as  discussed  In  the  appendix. 
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One  may  define  a  as  the  fraction  of  momentum  which 
is  used  to  excite  the  shock  wave  in  the  explosive  crystal* 
If  the  collision  is  perfectly  elastic 

e=  0, 

If  it  is  completely  inelastic 

«>a 

a  — - — __ 

mi+ms 

whore  mx  and  m2  are  masses  of  hitting  body  and  body  hitted, 
respectively..  In  our  case  a.  must  be  very  small  since  the 
collision  i3  nearly  elastic „ 

The  lattice  wave  propagates  with  the  velocity  of  sound 
ua  The  number  of  excited  oscillators, n,  Is  thus 

n  =  NAuAt 

where  N  Is  the  number  of  oscillators  In  the  unit  volume 
and  A  I3  t  he  area  vhich  is  hit  by  the  external  body. 

Prom  the  preceding  equations,  we  have 

At  -  ( aP/2NAuMwBx0)^ 

n  ~  ( aPNAu/2Ma)cx0 )  ^ 
if  At  «l/w,  and 

At  —  o  P  /KA'uMcox 
'  o 

n  =  aP/MwxQ 
if  At  p^l/u 
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For  example ,  if  the  hitting  body  has  the  same 
density  as  the  explosive  and  its  velocity  is  v 

aP/2FAuMx0  a  { v/u )  f. i/2xn ) 

^  l 

where  1  is  the  thickness  of  the  hitting  body.  If  v  is 
100  cm/sec,  A  Is  1  cm,  the  above  ratio  is  about  10^- 
sinca  u  Is  about  10 ^  crn/sec,  x0  is  about  10  “7  cm.  Thus 
if  c  »10“^  we  have  a  «At  case,  but  if  cr<S£l0"^-  it 
Is  a  «At  <T<  1  case. 

If  we  drop  a  mass  m  from  the  height  h  the  momentum 

P  is 

p  m(2gh)» 


thus 


n  =  (omNAu/2MwBx0)^2gh)1/^ 


if  wAt  «1,  and 


n  ~  { am/Mcox0 )  ( 2gh)  * 


if  wAt  »1. 

Now  the  energy  of  this  hitting  body  Is  mgh,  so  that 

W  «  Pmgh 

where  P  gives  the  fraction  of  the  energy  absorbed  by  the 
shock  wave.  It  is  expected  that  P  Is  about  0*5=  Using 
the  above  equations,  one  can  obtain 

*o  ■  «>* 


where 


I 
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if  (oAt<^l,  If,  on  the  other  hand  coAt>J>l, 


r  =  (P/aw)(g/2)^ 

if  coAt  »  1. 

We  see  easily  that 

At  =  2“^^2{cm/PMAuM®)^-\  If  wAt<<<l 
At  =  2osm/3MHAup '  if  o>At 


Suppose  h  is  1C  cm,  we  have  v  -  100  cm/see  which  gives 

the  previous  example ,  Since  w  is  about  10 -*2  sec“^,  a 

should  be  about  10“^°^  in  order  that  x  is  lO-"?  cm. 

o 

Prom  the  last  equation  above,  we  have  At  of  about 
lO"”10  sec.  Thus  we  have  a  consistent  picture  of  wAt»l 
case  here. 

Using  x0  «  Yh*  we  obtain  the  reaction  rate  II  for 
three  cases  considered  In  the  preceding  chapter. 


_  S*  k0 

2<JI  (h-h0): 


=  0 


bo 


if  h  >ho 
if  h  <h 

o 

if  h  >  h 
^  o 

if  h  <  h 

o 

if  h  >  h 

If  h  <  hQ 


wher  e 


h.  -  {a yr)B 

u  u 

In  Figure  A  the  corresponding  curves  are  shown* 

Since  the  velocity  of  sound  is  independent  of  the 
amplitude  the  above  quantity  essentially  gives  the 
number  of  reactions  per  unit  volume „  Implosion  will 
take  place  if  the  reaction  per  unit  volume  is  above 
certain  value  so  that  the  temperature  of  the  explosive 
can  be  abovs  the  explosion  temperature 0 

The  number  of  reactions  which  take  place  in  unit 
volume  is 

yMAt 

From  our  equations  for  At  we  see  that  At  is  independent 
of  h  in  any  case*  the  height  dependence  of  the  above 
number  of  reactions  per  unit  volume  is  given  by  the  K  part, 
namely,  by  our  equations  for  Kj. ,  Ks ,  Ka  for  three  typical 
cases,  respective  ly. 

The  explosion  can  oc  lur  if  the  temperature  reaches  to 
the  explosion  temperature  Te  by  these  reactions „  Assuming 
that  we  can  neglect  the  heat  conduction  during  At  the  condi¬ 
tion  of  explosion  is 

|  NKAt  ^  T0-T  , 

where  9  is  the  heat  produced  by  each  reaction,  C  is  the 
specific  heat,  T  is  the  original  temperature. 

For  a  given  size  of  the  explosive  the  explosion 
temperature  Te  is  given  as  discussed  in  the  preceding 
chapter . 


TO 


The  above  formula  gives  the  minimum  original  temperature 
T  necessary  for  explosion  as  a  function  of  the  height  h 
if  we  put  out  various  expressions  for  K  and  for  At  into 
it.  The  result  can  be  seen  in  Figure  B„  These  three 
cases  can  be  easily  distinguished  from  each  other  in 
this  case,,  Such  experiment  will  thus  be  a  good  test  to 
see  which  assumption  is  correct* 

Conclusions  and  Recommendations 

The  reaction  rate  theory  developed  in  this  text 
may  be  applied  to  the  description  of  explosion  initiation 
by  a  falling  body. 

A  preliminary  and  simplified  treatment  given  herein 
should  be  amplified  so  as  to  include  further  detail. 
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Reactiop 

Rate 


Figure  A„  Reaction  rate  as  a  function 
of  height  for  cases  1,  2,  3. 


Figure  B.  Minimum  original  temperature 
necessary  for  explosion  as  a 
function  of  height.  Curves  1, 
2,  3  correspond  to  our  cases 
1,  3,  and  3,  respectively. 
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Appendix 


CALCULATION  OP  THE  FORCE  DUE  TO  SHOCK  WAVE 

We  suppose  a  linear  chain  made  by  single  kind  of 
atoms  combined  together  with  force  constant  mo>8  •  If  we 
hit  one  end  of  this  chain  the  vibration  of  this  atomic 
system,  which  is  nothing  but  a  longitudinal  sound  wave 
or  lattice  wave,  will  be  excited  at  that  end,  The  wave 
will  propagate  with  the  sound  velocity  u  and  the  Inter¬ 
atomic  distance  r  at  distance  y  measured  from  the  end 
will  experience  a  displacement  x  from  its  equilibrium 
value  rc  as 

x  -  xQ  cos  -  t) 

If  the  excitation  mechanism 
is  such  that  it  keeps  pro¬ 
ducing  such  lattice  wave  at 
the  end  in  the  same  way 
throughout  a  period  At,  we  will  have  a  region  from  the 
end  surface  to  depth  y  =  ut  in  which  all  atoms  are 
excited  and  Interatomic  distances  given  by  the  above 
formula,  but  beyond  that  region  no  excitation. 

X  =  xQ  cos  w{^  -  t)  for  y  £  ut 
x  -  0  for  y  >  ut 


V=  'Co-tX 

CX-  .  — SO 


O.'toiryx 

O 


O 


t 


Figure  1,  Linear  chain 
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The  total  displacement,  namely  the  change  of  the 
length  of  this  chain,  at  time  t  smaller  than  At  will 
be  given  by  integrating  the  following  expression.  Using 
notation  N  for  the  number  of  atoms  in  unit  volume  and 
n  for  the  total  number  of  excited  atoms  after  the  time 
interval  At,  we  obtain  the  following  result,  since 

n  -  NAuAt 

The  total  displacement  at  fc  in 

x  “  NA  j  x  dy  -  NA  xn  /  eosw(X  -  t)dy 

■i  j 

=  nx„  Sinwt 
0  wAt 

The  total  force  due  to  the  displacement  is 

P  =  nMw8x0 

The  above  total  displacement  x  can  be  both  positive  and 
negative.  If  the  impact  is  so  strong  that  the  surface  cannot 
be  pushed  back,  the  negative  value  of  x  cannot  exist.  We 
assume  that  each  time  x  tries  to  be  negative  a  new  wave  is 
generated  at  the  surface  so  that  x  always  stays  positive. 

One  such  wave  is 

x  =  nxc  / 

In  this  case  we  have 

P  =  nM»Bx0 

?k 


for  example 


The  average  total  force  over  the  duration  time  At  is 

,  At 

<f>av  -  ^  dt 


At 


If  coAt  «  1  we  have  for  both  cases, .that  is  for 
F^sin  wt  and  (sin  cut  j 

<F>\V  = 

while  if  coAt»l  we  have 

<p>  =  nMw® xn  IpLVL^P,  ~  sMgx, 

'  AV  °  (aiAt  )E 

for  the  sin  fait  case  and 


ia 
(falAt) B 


/  p'>  £  nKcoExc  ^nMcogx0 

'  ''  AY  “  ti  aiAt  uAt 


in  the  case  of  |  sin  cot  |  „ 

If  we  take  P  for  the  sin.  cot  case ,  we  obtain 


P  =  NAuMx0 

which  means  x0  Is  about  10  ^  cm  under  ordinary  experimental 
setup.  Obviously  xQ  must  be  about  10 ”7  cm  in  our  picture. 
Therefor®,  we  do  not  tako  (f)  corresponding  to  Pm  sin  cot 
in  this  article  but  take  <^F>  corresponding  to  P^sin  cot 
instead. 
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ChaDter  IV 


THEORY  OF  DIELECTRIC  CONSTANT  AND 
LOSS  TANGENT  FOR  A  MEDIUM  CONTAINING  IMPURITIES 

Donald  G.  Burk  hard 

In  the  final  report  for  Contract  DA-44-009-ENG-3628,  1959 
we  reviewed  various  atomic  and  molecular  mechanisms  which 
give  rise  to  a  frequency  dependent  dielectric  constant 
and  loss  tangent. 

In  the  same  report  we  calculated  the  frequency 
dependent  dielectric  constant  and  loss  tangent  for  a 
medium  containing  impurity  centers.  The  principal 
motivation  for  investigating  this  situation  lies  in  the 
fact  that  newly  formed  ’’pure”  lead  azide  decomposes  with 
time.  Therefore,  if  dielectric  and  l.oa3  tangent  measure¬ 
ments  are  made  on  a  specimen  at  different  times,  one  is 
not  necessarily  dealing  tvlbh  a  pure  homogeneous  specimen. 

When  attempting  to  interpret  ralcrotfAve  measurements,  it 
is  therefore  of  importance  to  know  what  type  of  frequency 
dependent  behavior  to  expect  as  a  result  of  the  presence 
of  impurities.  The  measured  spectrum  will,  in  addition 
to  the  impurity  effect,  have  superimposed  on  it  the 
characteristics  of  the  molecular  polarization  associated 
with  the  pure  medium  and  with  the  impurity  centers. 

These  effects  may  be  incorporated  into  our  general 
formulas  for  particular  cases. 

As  we  shall  see  In  the  next  section,  the  exact 
calculation  of  the  dielectric  properties  of  an  Impure 
medium  is  analytically  rather  complicated.  One  may 
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substantially  simplify  the  calculation,  however,  by  using 
Lorentz’s  procedure  for  calculating  the  dielectric  prop¬ 
erties  of  a  medium  containing  polarizable  atoms  and/or 
molecules.  The  procedure,  as  is  well  known,  is  t o  remove 
an  imaginary  sphere  from  the  medium.  The  field  at  the 
center  of  the  sphere  and  associated  with  the  matter  outside 
of  the  sphere  is  then  computed  by  calculating  the  field 
associated  with  the  polarization  of  the  sphere.  To 
obtain  the  total  effective  field  at  the  center  of  the 
sphere  one  then  replaces  the  matter  in  the  sphere  and 
calculates  the  local  field  associated  with  the  individual 
charges  of  the  matter  in  the  sphere.  Knowing  the  total 
effective  field  at  an  atom  site,  one  may  then  express 
the  atomic  polarizability  in  terms  of  the  external  field 
and  finally  obtain  an  expression  for  the  complex  dielectric 
constant  of  the  medium. 

In  a  similar  manner,  the  Lorentz  sphere  procedure 
was  used  to  calculate  the  effective  dielectric  properties 
of  a  medium  containing  a  cubical  array  of  Impurity 
centers.  A  requirement  for  accuracy  of  the  calculation, 
however,  is  that  the  separation  of  impurity  spheres  be 
appreciably  greater  than  their  diameters.  It  is  not 
clear  just  what  the  limits  of  accuracy  of  this  model  are. 

It  was  therefore  considered  desirable  to  attempt  to  solve 
the  problem  exactly  by  solving  Laplace’s  equation  for  a 


7  7 


homogeneous  medium  containing  a  cubic  array  of  dielectric 
spheres „  It  was  founds  after  developing  the  above  des- 
cribed  Lorcntz-type  treatment  that  Rayleigh  had  calculated 
the  effective  conductivity  of  a  medium  containing  a  cubic 
array  of  conducting  spheres .  In  this  report  we  generalize 
Rayleigh^s  results  in  order  to  calculate  the  complex 
dielectric  behavior  of  such  a  medium,,  The  more  exact 
but  also  more  complicated  results  obtained  by  using 
Rayleigh  (is  procedure  are  then  compared  with  ours  in  order 
to  determine  the  range  of  applicability  of  the  latter „ 

That  is ,  we  compare  with  the  more  accurate  treatment  in 
order  to  determine  the  maximum  ratio  of  diameter  of 
impurity  to  separation  of  impurity  In  order  for  our 
simplified  procedure  to  yield  accurate  results. 

We  now  summarize  the  procedure  for  calculating  the 
complex  dielectric  constant  for  a  homogeneous  medium 
containing  a  cubic  array  of  dielectric  impurities.1  This 
will  be  done  first  by  following  Lorsntz»s  procedure 
and  will  then  be  carried  out  by  generalizing  Rayleigh«s 
method.  In  a  later  section,  the  effective  dielectric 
constant  and  loss  tangent  for  a  compound  medium  will  be 
shown  graphically  as  a  function  of  frequency  for  various 
percentages  of  "impurity"  constant  and  for  various  values 
of  the  dielectric  constant  and  loss  tangent  for  the 
medium  and  for  the  impurities. 
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Lorentg  Procedure 

Allow  the  supporting  medium  to  play  the  role  of 
free  space*  That  is,  the  supporting  medium  will  be 
described  by  the  dielectric  constant 


p.  0 

eo" 


e„K~ 

o  o 


*o<k; 


fib . 


toe 


E  cj 


We  now  treat  the  imbedded  spheres  as  ’’atoms”  In 
the  supporting  medium* 

The  Lorent z  sphere  is  shown  dotted  in  Figure  1*  If 

the  Imbedded  spheres  form  a  cubic  or  random  array  of 

dipoles  parallel  to  the  field  then  the  internal  field 

of  the  spheres  produces  no  effect „  In  this  picture 

we  now  want  to  correlate  e“  and  C  for  the  composite 

e  ® 

medium  with  the  macroscopic  properties  of  the  spheres 
and  the  supporting  medium*  The  atomic  polarization 
effects  usually  put  into  the  Olauslus-MosottI  equation 
now  are  those  associated  with  the  imbedded  spheres*  The 
throe  basic  equations  which  lead  to  the  Clausius-Mosotti 
equation  are  new 
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Figure  1,  Imbedded  spheres  K"  (K  ,  cr  )  in 
dielectric  K* ( ,  o!i0 ) ,,  Dgtted 
line  alluvia  Lorsntz  sphere  „ 


x*  b<  :.ng  fell-;:  effective  '%'A  the  overall  medium  is  then 
the  \i  o t  a  medium  consisting  of  the  imbedded  dielectric 
spheios  in  the  r,spaee"  of  the  supporting  medium,  o  is 
the  polarizability  of  the  individual  spheres  in  the 
supporting  medium »  The  above  equations  yield  the 
Clan  £  :.u  s  -  Mo  •  sotti  result 
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To  calculate  a  one  must  calculate  the  dipole  moment 
associated  v;ith  one  of  the  imbedded  spheres «  This  is 
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done  by  solving  Cor  the  potential  outside  of  a 
conducting  sphere  in  an  infinite  conducting  dielectric 
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First  term  contributes  a  potential  which  is  the  same 
as  that  of  c  dipole 
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] 0 1  /l|icOfc  ,  in  MS  units;  G  is  the  velocity  of  light  in 
free  space  ia  meters  per  second. 


*  Chapters  ‘Theory  of  the  Dielectric  Constant"  by  D.  £>. 
Burk-hard  In  Final  Report,  Contract  M-Wi-009"ENG-3628 , 
July  1959.  ' 
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t*/e  and  t^/tn  are  the  complex  dielectric  constants 
for  sphere  and  modinm.  The  coefficient  of  E  is  then  the 
polarizability  "ur<  ci  our  macroscopic  spheres  o 
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Putting  the  expressions  for  and  t:*  in  the 
equation  for  ?<  and  lotting 
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one  obtain.'. 


,  .  if  +  >1 K 

(J  «  .•*  ui .  imcn  :r4  >«.>■  >  -* 

Q”  +  ;jtoK 

if  we  now  pub  this  egression  for  c  In  the  Clausius- 
ios o t  b  1  h qu a  t i  on j,  one  has ; 


Oj 

”c. 


-■  j«Kf)  v  do  >)K  )  <r:t  ...B  _+  j(*K  Ci'/jK  "Bi 

c°  °  CT-Cl/^*  B  +  ju[K  ~(l/$C°B] 

»;here  B  na^  ~  nv,  v  is  the  volume  of  a  single 

Impurity  sphere  and  n  is  the  number  of  such  spheres  per 

unit  volume „  If  one  now  equates  real  and  imaginary 

>ar  ts  the  following  expressions  are  obtained  for  the 

effective  dielectric  constant ,  Kg r  and  conductivity  ^e. 
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of  r-bfi  composite  medium 


K 

e 


■a A  ;t*+«*k* }^I(a*<r‘>tt»KKv ;  b  £(a"?fc+wfcK  "£)be -t’k) 
3  9^  &o 

0“R  +  «®K*  -  ■?-  ajV  +  ojsKK 3 )  B-i-  -  (CT ' 8  +  jwaK''®)Bfc 
3  9 


■  e 

f'c 


“  s  [  (<T;^  )  •^•^•<r*+(0sKK  ' )  B  .~fer  *^w*K  '**  )  B* ;  ~«8K  (JCTK  *  <K)B 

9  o 


o  3 

<C  *  r  oo^K8 


CcTOr  +  u8KK?  )B  +  .|(or,ts  +  a>£KV!E)B® 


In  the  above  equations  one  now  has  K  and  CT"  expressed 
I  n  terms  of  F  0 „  'X0  ,  K., ,  <5"a  .,  B  v  on; 


K0  (K0,°o  „  Ka  ,  ,  B,  co) 

e  (Ko^ro>  Ks»^s»  B,  ») 

J i ayl el gfa  Pr o c edgy a 

Following  is  n  generalization  of  Rayleigh  's  procedure  ,1 
In  principle,  there  Is  no  restriction  on  the  size  of  the 
impurity  spheres.  The  results,  however,  are  expressed 
in  terms  of  a  series  expansion  so  that  one  must  be  careful 
to  carry  a  sufficient  number  of  terms  for  adequate 
i onvergenee , 

Philosophical  Magazine,  Vol„  33-3^  (1892)  pp„  lj.81-502. 


Consider  the  medium  as  divided  into  an  array  of 


rectangular  units  of  dimensions  a,  P,  and  Y  in  the  x,  y, 
and  2  directions,  respectively.  See  Figure# 


cylinder  way  be  expanded  in  the  aeries 


outside  V  “  AQ  +  (Axr  +  B1r“s)Y1  +  •  • »  +  (Anrn  +  Bnr“n‘'^)Yn 

inside  V9  «  0o  +  C^r  +  •••  +  CYrn+  ... 

0  A  n  n 


Yn  denotes  the  spherical  surface  harmonic  of  order  n. 
Boundary  conditions  on  the  surface  of  the  sphere  where 
r  -  a  ares 


VdV '  s;  dV 
dn  dn 


Y  denotes  the  ratio  of  the  complex  dielectric  constant 

of  the  sphere  to  the  complex  dielectric  constant  of  the 

medium.,  ... 

K 

T>  « 

K* 

o 

To  apply  boundary  conditions,  one  has 


dV  a  f  A 
dn  Ul 


2Bx 

r3 


)Ya  +  » •»  +  nAnrn  ^  -  i.£~3Jt  3..Y_  + 
n  p+2  n  n 


3i!pL  a  ^GYi  +  ...  |/nC  Y  r"*1  + 
dn  1  n  n 


Equating  coefficients  of  the  Y’s  for  r  =  a 


Ai  -  ~  =  M3x 

I  a-5  ' 

/  • 


( n+1 ) B 

A - S 

n 


na 


- 

2n+l  n 
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Likewise  from  V  -  V1 


A3. 


Ci 


A 

n 


+  -Sa—  *  cn 

a2»+l 


Prom  which 


n 


a2n*l  A 


l+V+l,  n 


n 


We  must  now  consider  the  limitations  to  be  imposed 
upon  Yn„  In  general. 


gS5  ^ 

=  2  $  (H3  cos  sq>  +  Kgsin  sip),  «»■ 


where 


<9®=  slnse(cosn"se  -  iB-AUaa=3l  oosn-»-2e  +  ...) 

2(2n-'J.) 


0  being  supposed  to  he  measured  from  the  axis  of  x 
(parallel  to  0}  and  «p  from  the  plane  of  xz„  In  the 
present  application  symmetry  requires  that  s  should  be 
even,  and  that  Yn  (except  when  n  -  0)  should  be  reversed 
when  (tt-6)  is  written  for  6.  Hence  even  values  of  n  are 
to  be  excluded  altogether „  Further,  no  sines  of  sip  are 
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admissible.  Thus  we  may  take 


Yi  -•  cos9 ,  «•*••••• 

Ya  =  cos^S  -  ^  cos9  +  Hz  sins9  cos©  cos  2<p, 

Y-^  =  cos^9  •=  cos^Q  +  Jj!_  cos6 

5  9  2i 

+  Ls  sin* 9  (cos30  -  I  CO30)  COS  2<p 

3 

+  sln^0  COS0  COS  1|<P  «  «  •  »  « 

In  the  case  where  P  Y  symmetry  further  requires  that 


He  =  0 »  Lg  —  0 


One  may  now  apply  Greenes  theorem: 


(tlSX  ~  ydU)  da 
dn  dn 


0 


to  the  surface  of  the  region  between  the  rectangle  ABCD 
and  the  sphere  P,  Within  this  region  V  satisfies 
Laplace ?s  equation,  as  also  will  U,  if  we  take 


U  =  x  =  r  cos 9 

The  applied  field  Is  taken  as  parallel  to  AD,  BC, 
thus  sides  DO  and  AB  are  equipotential.  Over  the  sides 


iJ  'T  -3t7  J”tr 

BG,  AD,  -g--  both  vanish.  On  CD,  represents 

the  field  strength,. 


E  » 


-:s- 


(tf"  +  JmSqKq 


where  J  is  current  density  and 


jii 

we_ 


Therefore 


dV 

dn 


ds  over  CD,  AB  represents 


jPr 


a 


o  o 


o,  The  value  of  remainder  of  the  integral 


1 


-  ( V~  ds ,  over  the  same  lines  is  -PYVi  where  Vx  is 
dn 


the  fall  in  potential  corresponding  to  one  rectangle, 
as  between  CD  and  AB. 

On  the  spherical  part  of  the  contour  over  which  Greenes 

theorem  is  to  be  applied 

U  ~  a  cos9,  22  ss  «  d£  =s  -  cosO  • 
dn  dr 
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a  o  o 


V  =  A0  +  (Axa+Bxa~B)Yx  +  ( A3a3+Baa"^)Y8  +  « • » 

H  -  -  ~  |  ijsa  -  “  (Ax-2Bxa-3)Yl  -  (3Aaa8~liE3a-5)Y3 

Substituting  these  expressions  into  Green«s  theorem 
yields 


■J 


U”ds 

dn 


■Xvf- 


[  A  +( Axa+Bxa"s  )Yx+(  Aaa3+Baa-4)Y8+ » •  •  ] 


Inserting  the  appropriate  expansions  for  Yx ,  Ya,  °°u  gives 
\(U^I  -  VdU)ds  «  ^a3( AX“*2Bxa"3)  (  |  (eos9)  sin9  cos0  d9  dtp 

J  dn  dn  JoJo 

-a3(3AaaJB-4Bsa^)  (  cos3q  -  1  cosS+Hasin^QcosBcosBqi) 

°  °  5sy 

XsInG  cos9  d0  dtp  -**•«  +  aBA0  (  i  sin0  cosB  d0  dtp 

fP  Pf  <t>4a 

+a“  ( Axa+Bxa“fl )  I  |  (cos0)  sln8  cos9  d9  dtp 
+aE(Aaa3+Baa  X  f  f  (cos-^0  -  —  cos9+H8sinB9  cos0  cos2<p) 

4>  Jo  $ 

X  sin8  cos9  d©  d<p  + 
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Simplifying , 


j(ud^  '  VH5ds  =  “(-Aia3-2Bi)Ix“(3Aaa^-4Baa“s)(Ia-  |  Ij 

9  +affiA0I0+(Aaa3+B1)Ix  +( A3a^+B3a-a ) (Is-  1  I1+HE 

zy  y 


«  CO 


where 


I  * 

O 


sinQ  cos8  d0  dcp  ~  0 


Ti  -  jj  ^  £cos8)  sinG  cos 9  d0  dcp  =  — (2it)  - 

nr 

(3ine0  cos 0  cosHcp)  sin9  cos©  d0  dcp  -  0 
*y 

la  =  C  f  ( c o s 3 0 )  sine  cos©  d8  dcp  ~  2(2?t)  = 

°o  C  cr 


etc , 


Thus 

J(IfS  =  V^)ds  =  |  ,  4|)(3Aaa^-4Bs 


a 


•••+(Aa.a3+B1)i2Jr»-(Jl2L  -  I  .  iffi)(3Aaa^-4B3a~®)  + 
3  5  5  3 


=  If-nBi 


Although  only  terms  up  to  Ya  were  used  above,  it  is  seen 

that  integrals  of  the  type  s  { oos 0)Yn( sine  d0  dcp)  will 

vanish  for  n  =  3>5»7»*°“*  due  to  the  particular 

numerical  coefficients  on  powers  of  cose  in  Yn  and 

also  to  terms  which  contain  cos  kcp,  k  =  2,l±,6'° •  * 
rant 

since  J  cos  kcp  dcp  =  0,  k  =  2,4>6’ * “ ’  . 


+H,Ia) 

Ia)+»»« 


) 


91 


Thus  fro®  8reen*s  theorem  one  has  finally , 


—  -  PrVi  +  li-ltBi  -  0 
J«£oKo 

The  potential  V  at  any  point  may  be  regarded  as  due 
to  external  sources  at  infinity  (by  which  the  flow  is 
caused)  and  to  multiple  sources  situated  on  the  axes  of 
the  cylinders o  The  first  part  may  be  denoted  by-Exo 
In  considering  the  second  it  will  conduce  to  clearness 
if  we  imagine  the  (infinite)  region  occupied  by  the 
cylinders  to  have  a  rectangular  boundary  parallel  to 
er  and  Even  then  the  manner  in  which  the  infinite 
system  of  sources  is  to  be  taken  into  account  will  depend 
upon  the  shape  of  the  rectangle,.  The  simplest  case, 
which  suffices  for  our  purpose,  is  when  we  suppose  the 
rectangular  boundary  to  be  extended  infinitely  more 
parallel  to  a  than  to  P.  It  is  then  evident  that  the 
periodic  difference  V*  may  be  reckoned  as  due  entirely 
to-Ex,  and  equated  to-Ec.  For  the  difference  due  to  the 
sources  upon  the  axes  will  be  equivalent  to  the  addition 
of  one  extra  column  at  +  x,  and  the  removal  of  one  at 
-  ,  and  in  the  case  supposed  such  a  transference  is 

immaterial »  Thus 

Vx  =*■  Eo 


*  It  would  be  otherwise  If  the  infinite  rectangle  were 
supposed  to  be  of  another  shape,  e*g,  to  be  square. 
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Placing  Vj.  —  Eo  in  the  above  result  obtained  by  use  of 
Greenes  theorem,  one  obtains 


J 

E 


4*Bi. 

«Pye 


»  £1  -  J*eJ£ 


0  o 


E  is  the  aoplied  field  so  that  ~  is  just  the  complex 

E 

dielectric  constant  of  the  medium,,  For  o>  -  0  it  is 
just  the  specific  conductivity  in  a  direction  parallel 
to  a  of  the  compound  medium,.  It  now  remains  to  connect 


E  with  Bl0 

Following  Rayleigh  we  now  calculate  Bx/E  approximately, 
limiting  the  treatment,  for  the  sake  of  simplicity  to 
the  case  of  cubic  order,  that  is  where  a  «  P  =  Y„ 

One  may  now  obtain  a  relationship  between  Bx  and  E 
by  equating  two  forms  of  the  potential  for  a  point  x,  y 
near  P  in  appropriate  limiting  oases.  For  example,  the 
potential  at  P  due  to  multiple  sources  Q  and  source 
-Ex  at  infinity  i3 


T  =  2T  Z  ( A  r  1  n  +  B  r'n“bx"  -  Ex. 
n  n  " 


Here  the  r }  and  the  arguments  of  the  spherical  harmonics 
are  referred  to  coordinate  system  located  at -each  Q. 

5 1  denote  the  origin  of  the  Qss,  If  the  Q,  sources 

are  allowed  to  recede  to  a  great  distance  only  the  terms 
in  B  will  contribute.  The  above  expression  does  not  include 
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the  effect  of  P  itself.  On  the  other  hand,  if  the 
sphere  f  is  excluded  (if  its  radius  goes  to  aero)  the 
potential  at  xy  is  given  by  our  original  expression 
with  the  terms  in  B  absent.  Thus  another  for®  for  the 
potential  is 

An  +  AirYi  +  A3r^V5  +  •«« 

Equating  the  two  forms,  one  has 

A0  +  Axx  +  Aa(x3  ~  3.  xr*)  +  • ° »  , 


in  which 

x’  s  x-f,  j1  «  y-^»  a9  “  z-^5 

where  £  are  the  coordinates  of  Q  referred  to  P, 

Hence  hy  differentiating  the  left,  side  once  with  respect 
to  x  then  setting  x=y  =  z  —  0,  Ax  -  His  obtained  and 
is  equated  to  the  corresponding  expression  on  the  right 
side.  In  similar  fashion,  differentiating  the  left  side 
thrice  with  respect  to  x  then  setting  x  -  y  -  z  ~  0  yields 
3!a3,  which  is  equated  to  the  corresponding  expression 
on  the  right.  Thus  *T 


91; 


x=y=2SB0 


At  the  origin,  x  -  y 
Ai  **  H  -  Bj.2  -  "gj. 


Since  x5  =••  -J,  y5  -  ,  z*  ”  at  the  origin,  so 

also  is  r°E  "JE  +  *1 8  +'  f3*  and  thus 


4.  /i l  u 
d3  ( p)  ft* 


If  we  let  ,/¥ 


J/X  then 


37*-/** 

2/P 

t 


-  2_pB(/a 


where  PE  (.#)  Is  the  second  order  Legendre  polynomial . 


/°7l33“-  jy^) 


Z5 


ill 


35^f  -  3o£.  * 
^>U  ^ 


3 


-  JL  /  ) 

*4  8  J 


Also 
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c  %  tl53j^7  -  ?0^3/J"9  +  63J^’11] 

*  -  3/r7(23l3L  ~  315  *  105  It  -  5) 

f>\  {* 
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Hence 


A3  -  -UBa.  -  8Ba  2/5“7p6(^  -  ••• 

A  similar  procedure  will  give  Ag,  Ay ,  etc.  Each  of 
the  above  summations  is  over  all  the  other  spheres  Q, 
I.e„,  over  all  coordinates  y »  /  of  the  formic,, 
me,  no  where  /f,  m,  n  are  positive  and  negative  integers 
except  0,  0,  0.  Let  these  sums  be  denoted  by  Se ,  S^, 
etc.  so  that 

2/~n°'X  P„U0  *  o"”"1  SR 

Since  Ar  »  ®n ~2n+l  due  ^oundary 

conditions,  we  have  the  followings 


Aa  -  Ba  p  .liB^c-Ss.  )  -8B.(c"V) 


and  thus 


Ai  -  H  “ 


Bi  (limi  ..  H 

(l-^)a^ 


*2BJL(e'“38.)  -  |  B3(o“5Sk) 


or 


H 


*^*5  +  a*’-3*.*  !  «-S 


.l4Bl0-5S.  .  <i^l£ 

4  (T^T 


if  we  neglect  the  last  term  in  the  A3  expansion 


Therefore 


Finally,  the  specific  conductivity  is 

(where  for  simplicity) 


IprtBi  ■ 

1  -  -  »  1 


Una3 


oPYH 


;3(tf 


3 


J4.1t  a  3/^,3 


»  1 


f±2  +  2Ib\\  .  4  /if °/W)  sf 

1-2/  [a/  8  5  (0/  li¬ 


lt  remains  to  calculate  the  sums 
sum  Sa,  we  have 


For  the  first 


S2  =  2^“3Pa(/,)  * 


It  has  been  assumed  that  the  extension  of  the  medium 

is  infinitely  more  parallel  to  the  c  direction  than 

either  P  or  T  so  the  space  to  be  summed  over  is  a 

rectangle  bounded  by^5*  -  00 ,  ^ s  tV)  Js  where  v 

will  ultimately  increase  without  limit. 

Consider  first  the  space  bounded  byjf'=^"^~  ~v, 

a  cube  centered  on  the  origin.  Since  the  spheres  Q 

are  eubically  symmetric,  2  2L.  2  2—  =  2  so  that 

p$  p$ 

Sg  s  2  — p-  ~-s—  s  0  over  this  space,, 


Consider  next  the  space  bounded  by 
t  cxp  ~va  We  assume  that  in  this  space, 

f is  sufficiently  large  so  that  the  summation  can  be 
replaced  by  an  integral „  Hence 


Sg  =  -  A  2 

8  dJU? 


1 

2 


r 


y  j 


where  is  integrated  from  -  co  to  ~v  and  +v  to  +  ao  „ 

Hence  :|sa  will  be  obtained  by  integrating  from  +v  to  +  oo. 


-V 


»  2v* 


A 


(**+}*)(  2va-»5E)1/2 
Let  ^  ~  v  tan9,  then  d  j}  51  v  see80d9,  v®+^£  »  v*  see* 

(2v*+J8)l/2  *  v(2+tana0)l/2, ^  0  ~  J 


^  a  _v  — ? 0  s  ~jj-  so  that  the  above  integral  becomes 

■  (  2d9 

J  /  2+ tan* 9 


Sg  * 


9 


Let  x  ~  sinG,  then  dx  £“  eosQ  d8,  tan^Q 


.x2 

i=xs 


Finally  if  we  let  p  25  and  neglect  the  term 

(a \10  a3 

containing  /— j  as  a  first  approximation  to  the 
specific  conductivity,  -we  obtain 


In  order  to  carry  on  the  approximation  we  must 
calculate  S^,  etc.  An  approximate  value  of  may 
be  calculated  by  direct  summation  from  the  formula 


% 


3^|-3oT 


We  may  limit  ourselves  to  the  consideration  of  positive 

’  i  • 

,  are  finite  is  repeated  in  each  octant,  that  is 

8  times.  If  ono  of  the  three  coordinates  vanish,  the 
repetition  is  fourfold,  arid  If  two  vanish,  twofold. 


and  zero  values  of  ^ 


J"  .  Every  term  for  which 
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The  following  table  contains  the  result  for  all 
points  which  lie  within/*®  -  18-  This  repetition  in  the 
case,  for  example,  of  f>*  =  9  represents  two  kinds  of 
composition.  In  the  first 

p  ~  28  +  2b  +  Is  K  9  s 

and  in  the  second 

f>*  s  38  +  o2  +  0+  ~  9. 

The  success  of  the  approximation  is  favored  by  the  fact 
that  P  vanishes  when  integrated  over  the  complete 
sphere,  so  that  the  sum  required  is  only  a  kind  of 
residue  depending  upon  the  discontinuity  of  the 
summation. 

The  result  is 

Si.  -  3  ......  o  .. 


- , - 

1 

J _ 

IQQI 

1 

0,0,3 

9 

+  ,0lkk 

2 

-  .  309lp 

0,1,3 

10 

+.02k3 

| 

3 

-  .1996 

1,1,3 

11 

+  .0075 

0,0,2 

k 

+  ,109k 

2,2,2 

12 

-.0062 

0,1,2 

5 

+  .3501 

0,2,3 

13 

-  .0015 

1,1,2 

6 

-  .0397 

1,2,3 

lk 

-.0095 

0,2,2 

8 

0,0,  k 

16 

+  .00  3k 

1,2,2 

9 

-  .0277 

2,2,3 

,  17 

-  .0061 

!p,i,^.\ 

17 

+  .0085 
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As  indicated  earlier  the  preceding  several  pages  lor 
evaluating  E/Bx  follow  Rayleigh.  Details  have  been  added. 

Referring  to  our  formula  for  complex  dielectric 
constant  (talcing  a  —  £  =  Y) 


K*  -  [1 - kE — ] 

6  o3  E/Ba 


One  may  now  insert  in  the  above  the  following  expression 
for  E/Bi 


l 


lit:  a t3 

3  tr3 


(l-e*/ef)(3.11)g 

?z*/e*+l-24/34 


We  now  wish  to  compare  the  above  expression  for  K0 

with  our  form  for  K'*,  that  is  with  the  Lorentz  method 

for  obtaining  K  .  In  the  latter  case  we  have  (in  c  g  s) 
© 

units 


1  +  _  l+f|3tn0o  2/3 

l-lf-nnoo9/ 3  l-4mn0o  s/3 


In  this  case  cr 1  is  the  polarizability  pf  an  impurity 

t 

sphere. 

y 

First  note  that  if  Ipnn0c  /3  is  much  less  than  unity 
one  can  write  the  above 

K*  ~  1  +  ipin  c" 
e  ^  o 


102 


Also  note  that  this  is  the  form  one  would  obtain 

for  K"”'  if  it  is  assumed  that  the  local  field  is  equal 

to  the  external  field.*  That  is,  it  is  the  form  of  K* 

e 

when  one  neglects  the  effect  of  the  spheres  on  each 
other . 


o8  = 


1 


i Sv  * 


i  +  2 1  ;,*/£- 

0  s 


Putting  this  in  our  approximate  expression  for  K  ,  one 

e 

has 


* 


K  -  1  +  4nn0a: 
0 


„  «•  .  * 

1  "  Vh 

1  +  2B**/e* 

1  0  3 


If  one  approximates  E/Bi  by 


r 

1 

+  1 

a3 

eo/es  “  1 

-  -*•* 

P  a  XE 

Also  P  »  n„eE  =  n  trE 

o  eff  o 

if  local  field  is  equal  to  applied  field  E. 

Therefore 

X  =  n  a* 

o 

and 

K  -~  1  +  IjiiX  -  1  +  lj.nn  a5 
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In  the  Rayleigh  general  result 


one  has 


[1  -  ..iff— -J Kg 
o3e/Bj. 

1  +  1  ~  En/es 

1  +  aeg/eg 


1  _ 


=  n0  ~  number  of  spheres  per  unit  volume.  Thus 


employing  the  first  term  only  in  the  Rayleigh  expansion 
is  equivalent  to  theLorentz  treatment  when  the  local 
electric  field  is  equal  to  tho  external  applied  field. 
Since  this  is  not  true  in  our  case,  the  Rayleigh  result 
using  only  the  first  term  in  the  E/Bx  expansion  results 
in  a  poor  approximation. 

Now  examine  the  effect  of  including  the  second  term 
in  E/Bx.  The  first  term  is  simply  the  reciprocal  of  the 
polarizability,  o\  of  a  single  impurity  sphere.  Thus 

e/bx  «  [  -  JL  +  tea] 

o5  3 

Putting  this  in  the  Rayleigh  expression  for  K 


K*  -  1 
e 


1_  +  ^no 


=  1  + 


Jj.xn  o' 


This  is  our  general  result  using  the  Lorentz  procedure. 
Thus  one  can  say  that  the  inclusion  of  the  first  two 
terms  in  the  expansion  for  E/Bx  in  the  Rayleigh  method 
is  equivalent  to  the  Lorentz  procedure.  The  third  term 
in  Rayleigh') s  procedure  will  represent  the  improvement 
to  be  gained  from  the  more  exact  procedure „  In  view 
of  th8  above  results,  it  is  now  convenient  to  write 
E/Bj  as  follows 


a^(3,ll)c 


2*o/«S 


3(l*2e;/e*>]  J 

_  +  _  4-^PqC 

cr?  3  -  3 

lpnoC 

-  stands  for  the  third  order  correction  terms » 


E/Bi  = 


f 


o 


Then  the  Rayleigh  K*  may  be  written 

e 


K*  =  1  - 
e 


^no 


•j  ^itn0 
“  j  +  —  (1  -  o) 

a  3 


1  + 


c)a‘ 


Thus  tbs  Rayleigh  treatment  is  equivalent  to  writing  for 
the  effective  fisld  (in  cgs ) 


E 


eff 


=  E0  +  YP 
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where 


Y  s  krc(l  i  e] 
3 

rather  than  simply  Y  =  itE, 

3 


Since  there  is  no  advantage  in  displaying  the 
polarizability  </  in  c,  we  write  e  as  follows  for 
computational  purposes , 


c 


^(3„ll)e 


Cl+l^/3e" ) 


,10 


.10 


_1 _ 

n 

o 


a]£  Jsl)10'3  m  (n  a3,W/3 

10  3/  '  oa  ‘ 


B 


10/3 


-i  3n0v 

where  B  =  n0 aJ  -•  — j-— -  .  v  is  the  volume  of  a  single 
impurity  sphere  and  n0,  as  before,  the  number  of  such 
spheres  per  unit  volume. 


c  =  6,190 


/Ve 

■eo/es 


B?/3  »  2.168 


(n0v) 


7/3 


Prom  this  result  it  is  clear  that  the  accuracy  of  the 
Lorentz-type  treatment  depends  not  only  on  the  fractional 
volume  occupied  by  impurity  spheres  but  it  also  depends 
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on  the  ratio  of  the  medium  to  impurity  dielectric 
properties  as  was  to  be  expected. 

We  examine  c  for  the  limiting  cases  of  a  pure 
dielectric  impurity  (CT -  0)  and  for  a  pure  conducting 
impurity. 

Dielectric  case 


c  =  61.90 


1~K 


o' 


1 


l+ilK0/3Kc 


B 


7/3 


2.188 


i+te0/3Ks 


Conducting  spheres 


61.90 


l+I|(T0/3cr' 


J/3 


2.188 


(n0v) 


7/3 


Intermediate  cases  must  be  examined  by  equating 
real  and  imaginary  parts  of  our  general  expression  for 
K*.  Results  will  also  be  frequency  dependent. 
Expressions  for  KQ  and  <7^  must  then  be  compared  with 
those  obtained  from  the  Lorentz-type  treatment. 

We  shall  use  the  criterion  that  o  must  be  less 
than  0.1.  Figures  3  and  I4.  show  K0/Ks  as  a  function  of 
the  fractional  volume  n0v  for  c  *  -0.1  and  c  =  +0.1 
respectively.  Note  that  a  large  variation  of  K0/Kfl 
(2  to  greater  than  10)  is  permitted  when  the  impurity 
content  is  as  great  as  2/5  of  the  volume  of  the  sample . 
Similarly,  when  K  /K  is  less  than  0 .^.fractional  volume 

OS  / 


if, 


Figure 


Kc/Ka  as  a  function  of  fractional  volume 
nQv  for  c  =  +0.1. 
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of  the  allowable  Impurity  may  range  as  high  as  20  to  lj.0 
percent  for  the  Lorentz  model  to  remain  valid.  Thus  one 
can  conclude  that  for  reasonable  ratios  of  the  dielectric 
properties  of  the  sustaining  medium  to  the  impurity 
medium,  the  Lorentz  model  may  be  applied  with  considerable 
accuracy  even  though  the  fractional  volume  of  impurity 
becomes  rather  large  ,  that  is,  of  the  order  of  20  to  I4.O 
percent . 

Conclusions  and  Recommendations 

The  dielectric  properties  of  a  medium  are  sensitive 
in  the  microwave  rsgion  to  the  impurity  content.  If 
time  and  interest  permits,  it  would  be  worthwhile  to 
develop  the  typo  of  model  described  in  this  chapter  so 
as  to  include  radiative  effects  which  will  occur  at  the 
very  high  frequencies.  In  the  event  that  further  determin¬ 
ations  are  made  of  the  dielectric  constant  and  loss  tangent 
of  lead  azide  in  the  microwave  region,  one  should  look  for 
the  characteristic  behavior  displayed  by  the  curves  in 
this  chapter  in  order  to  determine  whether  impurity 
effects  may  be  significant. 


110 


(olT/iUV  w;j*v;  or?  ttT:!  DYRreoTp'i'g  rBONF.RT**"^  no  ,,y 
MEDIUM  >ifl£UK  CONTAINS  STHSRIOAL  DlfiJiEiGTfi XC  P/iiilM-J.uWo 

Do  !■’  *  Schwalw 

Consider  a  uniform  dielectric  mod*  urn  „  hereafter  designated 
as  the  support,  mod' am  and  denoted  by  a  aubserlpt  ,  which 
contains  a  cubical  or  random  array  of  Imbedded  spherical  particles, 
hereafter  doai  fixated  as  the  impurity  medium  and  denoted  by  the 
subscript  .  The  support  medium  will  be  characterized  by  a 
dielectric  constant  K0  and  a  conductivity  0O  while  the  impurity 
modi urn  will  be  characterized  by  a  dielectric  constant  Ks  and  a 
conductivity  CT.  „  rt  is  desired  to  obtain  the  resultant  behavior 
of  the  combined  media  which  can  bo  characterized  by  some  "effective” 
di  electric  constant  Ke  and  some  "offootivo”  conductivity 

In  a  paper  or.  the  THEORY  OE  THE  DIELECTRIC.!  CONSTANT >  D„  O. 


Burkhard  has  derived  an  extended  vorsi on  of  the  usual  Clausius- 
Moootti  equation.  tfe  summarize  that  derivation  here.  Three 
basic  equations  of  sloe fcromtvme tic  theory  lead  to  the  Clausius - 


Moootti  result,,  They  are 


P 


V 


EE 


eV( 


ri¬ 


ot.  Eett 
P 

36.'* 


t  1 y  v' 

whe^e  ~  ^0^0  is  tho  complex  fctivl  ty  of  the  support 

medium,  Xg  is  the  effective  complex  susceptibility,  oL  is 

the  polar! sabi  11. ty  of  a  spherical  impurity,  and  Ez  is  the  local 

field  in  the  vicinity  of  an  impurity  sphere.  The  EE  cff  is 

obtained  by  calculating  the  field  In  a  spherical  cavity  du©  to  the 

polarization  of  the  support  medium  and  the  field  due  to  other 

impurity  spheres  in  the  cavity.  Since  it  is  assumed  that  the 
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impurity  spheres  form  either  a  cubic  or  random  array,  this  latter 
fiold  vanishes. 

Elimination  of  "P  ,  P0  »  and  E:  eff  yields 
y  *  _  °l 

I  -  ^  ^ 

However,  to  complete  the  derivation  we  need  to  use  the  conventional 
relations  that  exist  between  the  ordinary  dielectric  constant  and 
conductivity  and  the  "complex  dielectric  constant,  namely 


x  *  -  C  -  1 

e  c 


Thus  we  obtain 

__  * 
€ 


I  ^ 
€c 


770c4 

!  -  ^o°V3 


+  jw  Kc 


*  jw  K 


i  +  3  fy>oi. 

I  -  i  A>«i 

\  f<+  jy>A 
/Li  - 


Now  It  remains  to  calculate  the  oi  of  one  of  the  imbedded 
spheres  to  complete  the  derivation. 


Assume  that  the  volume  loading  is  small  so  that  the  problem 
of  finding  reduces  approximately  to  the  problem  of 

determining  Os  for  a  single  dielectric  sphere  in  an  infinite 


medium.  In  this  case,  solution  of  Laplace’s  equation  yields  for 
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the  potential  outside  the  sphern 

/  c*  r'*  r-* 

I  CL 


[  ae'*  t  6/  r‘ 


Ee«  2 


which  consists  of  a  dipole  term  and  a  uniform  field  term.  The 

dlrsnlA  nntfint.l  nl  -1  a  C(j  —  !-*>£  OCOit} 


dloole  potential  is  Cp 

dipole 


so  that  the  lipole  moment  is  thus 

€,  -  2  f- 


and  so 


3€'  *  t  €  * 


€  *-  e#* 

co _ 

£€/'+  €* 


<ae„'%e,‘  m  eJ'/e/ 

where  CL  la  the  radius  of  an  impurity  sphere.  Substituting  this 
value  of  oC  Into  the  Olausius-Hosottl  equation  yields 


have  been  substituted. 
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TMs  express!  or  can  now  be  reduced  to  obtain  explicit 


equations  for  CT  and  which  is  the  final  desired 


result.  Tf  we  write 


B-  ria 


<r- 


c r'*  i 


(o;  f  3%)  ;  k 
1  (°s "  t 


then 


"(I 


'*  1  3 
CL 


(o;-o;)+ja>(k^Ko)  3 


_  cr  +  i  co  k  s 
-  - £ — _  a 

cr  t-  j  oj  k. 

so  that  the  Olausl.us-Mosotti  equation  becomes 

a. 

3 


1 ♦><. - (| *>«.)  it.iJvLkk^ 


a  Decomes 

\  cr+  jju  k  7 


i  »0a3 

V  <r  +  j  a>  *  / 


(To 


f  jcoJ 


v  2^1 

r  (T^JUIK  \ 

\  1  +  3  i5  i 

V<r  +  J“*  j  . 

'  l~  is 

n  jwi 


5  (r»f3r  )*jm(Ki-§3K'j  ,u||  (c~t  j3<r)+ju'(«  +  t3fcl) 

€-  C<r-i3ir')+jW(k-i3K')  J"  ”  (V-|3r0rr(k'-iBi<') 


Tf  subscript  o l  denotes  the  form  X  +•  f and  subscript 
j3  denotes  the  form  X"  {^X  ,  then 


°o  /  0^  T  ,j  CO 

€o\^  t  jcoV 


IjojKJ  --at 


<r*  4-  gj*k* 

li  P 


Separating  and  equating  real  and  imaginary  parts  gives 

r*  =  -g-C^O-.^XK^-coXC^k^-  ZK') 


<5*  *  tozK^ 


Kb.  f 

V  *  O)2^2* 

Keplacing  0^  by  CT  i  §3cr  ,  0"  by  CT- ~13cr’  »  by  R +  §  15  k 

by  R  -  and  separating  out  the  Ui  dependent  terms  we 

obtain  as  intermediate  results  v.  \ 

K.  (0-*  |-Br'Xir-  iTo-'j  +  -^(oV-ir'*)  t m  ~ ) 

ttT  *  (jx~  ^3(r')z  +  of-  (k-  i$K,)a‘ 


o;  __  (<r+-  |3<r')(<r-  l^r')“^~^tu^'-<TK)4-cJ^Kt"BK](K--3  jBk; 

(^r- ^o-')2,  +  t*ja  (k- i’BK7)1 

Now 

CT*-  §^<r'  =  ^  [(i+  §3)0;  +  (5-f3)o;J 

r-  i  'B<r'  -  1  [  (I  “  *•  (a  i-  je)  a; 

f"* 

(TR  (TK  ®  —  (TK  “  <r  k 

e«  L  41  5  *  UJ 

However,  33  =  7loa3  -  ^Jj&l  =  ^~T  **oU 

*  — ^ — {  fractional  volume  occupied  by  Impurities  ) 
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To  be  const stent  wi fch  the  model,  the  fractional  volume  of 
Impurities  must  certainly  not  exceed  25$.  Even  then,  -t>  —  O.OC, 
so  that  the  error  Introduced  into  CT+'^'Bcr  ,  CT- ^  3cr  ,  K+'|T3k<  , 
and  K -  ^'Bk'  by  neglecting  3  relative  to  unity  will  not 
exceed  4$.  For  plotting  purposes,  this  error  is  negligible. 


Thus  we  write 

Kc 

-  rs* 

Ko 


5, 
c r„ 


,  „  /  /  t  \ 

+  -i^C™  -*** ) 


<r 


crw  t 


■it-  *■ 


>6dV~ 


CT2'!  ^ZK 


or 


Kc 


,  ,  “£0; 
I  * 


e<,K'0 


" ‘1=-,';  *■  4<E *)  <•  e>%f H k. K, . Ik.1) 


I  +  em 

O 


cr0  L(^’*V "4%  +  e0za)X^  *■  4  k*  k 0 14/<  *  ) 


Now  introduce  the  following  ratios: 

<r£  .  0  _  Ks 


Kt 


Ko 

(Tc 


Then 


^  _  _ 33(&~#) 

K0  '  (7-^7'H)  -j-  +4) 

1  + 

*4)  +- 


a. 

cr0 


As  33  ®-nd  333  <&  43  to  the  same  approximation 

previously  made,  i.e.  neglect  of  3  compared  to  unity,  tae 
above  two  ©xoress-'  ons  will  reduce  to 


J4 

Ko 


I  i- 


3-3  e 


(<pzh4&  M)  +-  e*p*cj x(02h43  tT) 
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Rnd  .  3B 

c r0  '  +  {f-+44+4)+  €?f*?(<9z*<s*A~)" 

This  5.a  the  form  which  will  be  moat  convenient  for  plotting.  In 
the  accompanying  graphs,  ^£/\<0  have  been  plotted  for 

the  selected  values  of  0  =  10,  100,  1000,  and  10,000  and 

<J&  =  10,  10^,  10  A  10?,  and  10^  over  a  range  of  CO  from 
CO  =  1  to  CO  -  lO^/seo.  The  value  of  p  has  been  taken  to 
be  p  -  =  5  *  10-1*5  sec2. 

^  1>  ft"  I 

It  can  be  seen  that  — £  — *  I  as  oo  --*-ao  ,  - >-  /  as 

Ko  CTo 

CO — 0  ,  >.  _^LE£L__  as  C<J  — >*  o  ,  and 

K*  0*t4s*+4 

°e.  a  I  33#  as  o>  — *» oo  .  Thus  the  general  shape  of 

<r0  "**  9$*  +A<i> +-”4 

all  curves  will  be 


If  we  return  to  the  Intermediate  results,  we  can  evaluate 


the  limits  of 


V* 


and 


A 


as 

Thus 


oo 


which 


corresponds  to  metallic  impurities, 

Jfew.  K£  jllm,  K+|Bk  _  jfcw. 


Ki~>co 
CO  -V  oo 


co~*  o 


H  ~»cd 


k--;-bk' 


It  §3 


l+S 


e.  __  i 
Kt  ?  1?  k 


+  3 
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In  summary,  the  imbedded  sphere  model  with  sphere 
concentration  being  small  leads  to  a  frequency  dependent 
effective  conductivity  and  effective  dielectric  oonstant. 

In  addition,  it  shows  that  -  /  and  -S-  —  I  are 
directly  proportional  to  the  volume  of  the  impurity  con¬ 
tained  In  the  support  medium. 

Conclusions  and  Recommendations 

The  presence  of  impurities  in  a  dielectric  medium 
affect  the  frequency  dependents  of  the  dielectric  constant 
and  loss  tangent.  Curves  showing  this  effect  are  presented 
in  this  chapter  for  various  values  of  the  parameters 
characterizing  the  sustaining  medium  and  the  dielectric 
impurities. 
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